HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoFeminosAnalyseNoiseMonitor.cxx
Go to the documentation of this file.
1 //
2 // File HarpoFeminosAnalyseNoiseMonitor.cxx
3 //
13 #include "HarpoConfig.h"
14 //#include "HarpoDetSet.h"
15 #include "HarpoFeminosEvent.h"
16 
17 #include "TROOT.h"
18 #include "TStyle.h"
19 #include "TApplication.h"
20 #include "TFile.h"
21 #include "TCanvas.h"
22 #include "TLatex.h"
23 #include "TGraph.h"
24 #include "TF1.h"
25 #include "TH1F.h"
26 #include "TH2F.h"
27 #include "TLinearFitter.h"
28 #include "TMath.h"
29 
30 #include <cstdlib>
31 #include <cstring>
32 #include <cassert>
33 #include <iostream>
34 
36 
38 {
39  //int nd; // number of detectors
40  //HarpoEventHeader *hdr = anaEvt->GetHeader();
41 
42  // assert(hdr != NULL);
43  // hdr->print();
44 
45  // for (nd = 0; nd < gkNDetectors; nd++) {
46  // // if (anaEvt->isdataExist(nd)) {
47  // HarpoDccMap *plane = anaEvt->GetDccMap(nd);
48  // if (plane != NULL )plane->print();
49  // }
50 }
51 
53 {
54  nEvents++;
55  std::cout << " Processing Event " << nEvents<< std::endl;
56 
57 
58 
59 
60 
61  TCanvas *c = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("cHarpo");
62  if(!c){
63  c = new TCanvas("cHarpo","FEMINOS: All Events cumulated",0,0,1000,800);
64  c->Clear();
65  c->Divide(3,3);//,0.001,0.001);
66  }
67 
68  TCanvas *c2 = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("cHarpo2");
69  if(!c2){
70  c2 = new TCanvas("cHarpo2","FEMINOS: Current Event",0,0,1000,800);
71  c2->Clear();
72  c2->Divide(3,3);//,0.001,0.001);
73  }
74 
75  gStyle->SetOptStat(0);
76 
77  TH1F* h1Ev1Ch = new TH1F("h1Ev1Ch",";time [tb]",508,0,508);
78  TH1F* h1EvAllCh = new TH1F("h1EvAllCh",";time [tb]",508,0,508);
79  TH2F* h1Ev2D = new TH2F("h1Ev2D",";time [tb]; channel",508,0,508,304,0,304);
80 
81  TH1F* h1EvFFT1Ch = new TH1F("h1EvFFT1Ch",";freq []",508,0,508);
82  TH1F* h1EvFFTAllCh = new TH1F("h1EvFFTAllCh",";freq []",508,0,508);
83  TH2F* h1EvFFT2D = new TH2F("h1EvFFT2D",";freq []; channel",508,0,508,304,0,304);
84 
85  TH1F* h1EvPed1Ch = new TH1F("h1EvPed1Ch",";pedesteal [ADC]",500,0,500);
86  TH1F* h1EvPedAllCh = new TH1F("h1EvPedAllCh",";pedesteal [ADC]",500,0,500);
87  TH2F* h1EvPed2D = new TH2F("h1EvPed2D",";pedesteal [ADC]; channel",500,0,500,304,0,304);
88 
89  // Int_t npix = 0, nch = 0, ntb = 0;
90 
91  HarpoDccMap * m = anaEvt->GetMap();
92  std::cout << "Event " << nEvents << " " << m << std::endl;
93  if ( m != NULL ) {
94  // for(j=0;j<NCHAN;j++) // Channels
95  //Double_t qTot = 0;
96  TH1F* hTmp = new TH1F("hTmp","",508,0,508);
97  for(Int_t i=0;i<NADC;i++){ //Time bins
98  Double_t qTotT = 0;
99  for(Int_t j=0;j<NALL;j++){ // Channels
100  // Double_t q = m->map(j,i);
101  Double_t q = m->GetData(j,i);
102  h2D->SetBinContent(i+1,j+1,h2D->GetBinContent(i+1,j+1)+q);
103  h1Ev2D->SetBinContent(i+1,j+1,q);
104  if(j == 100){
105  h1Ch->SetBinContent(i+1,h1Ch->GetBinContent(i+1)+q);
106  h1Ev1Ch->SetBinContent(i+1,q);
107  hPed1Ch->Fill(q);
108  h1EvPed1Ch->Fill(q);
109  }
110  hPed2D->Fill(q,j);
111  h1EvPed2D->Fill(q,j);
112  qTotT += q;
113  }
114  hAllCh->SetBinContent(i+1,hAllCh->GetBinContent(i+1)+qTotT);
115  h1EvAllCh->SetBinContent(i+1,qTotT);
116  hPedAllCh->Fill(qTotT/288);
117  h1EvPedAllCh->Fill(qTotT/288);
118  hTmp->SetBinContent(i+1,qTotT);
119  }
120  TH1* hm1 = 0;
121  hm1 = hTmp->FFT(hm1, "MAG");
122  for(Int_t i=2;i<512;i++){ //Time bins
123  hFFTAllCh->SetBinContent(i+1,hFFTAllCh->GetBinContent(i+1)+hm1->GetBinContent(i+1));
124  h1EvFFTAllCh->SetBinContent(i+1,hm1->GetBinContent(i+1));
125  }
126  hTmp->Delete();
127  hm1->Delete();
128 
129  for(Int_t j=0;j<NALL;j++){ // Channels
130  TH1F* hTmp = new TH1F("hTmp","",508,0,508);
131  for(Int_t i=0;i<NADC;i++){ //Time bins
132  // Double_t q = m->map(j,i);
133  Double_t q = m->GetData(j,i);
134  hTmp->SetBinContent(i+1,q);
135  }
136 
137  TH1* hm = 0;
138  hm = hTmp->FFT(hm, "MAG");
139  for(Int_t i=2;i<512;i++){ //Time bins
140  hFFT2D->SetBinContent(i+1,j+1,hFFT2D->GetBinContent(i+1,j+1)+hm->GetBinContent(i+1));
141  h1EvFFT2D->SetBinContent(i+1,j+1,hm->GetBinContent(i+1));
142  if(j==100){
143  hFFT1Ch->SetBinContent(i+1,hFFT1Ch->GetBinContent(i+1,j+1)+hm->GetBinContent(i+1));
144  h1EvFFT1Ch->SetBinContent(i+1,hm->GetBinContent(i+1));
145  }
146  }
147  hTmp->Delete();
148  hm->Delete();
149 
150  }
151 
152  hFFT1Ch->GetXaxis()->SetRangeUser(-1,255);
153  hFFTAllCh->GetXaxis()->SetRangeUser(-1,255);
154  hFFT2D->GetXaxis()->SetRangeUser(-1,255);
155  h1EvFFT1Ch->GetXaxis()->SetRangeUser(-1,255);
156  h1EvFFTAllCh->GetXaxis()->SetRangeUser(-1,255);
157  h1EvFFT2D->GetXaxis()->SetRangeUser(-1,255);
158 
159 
160  c->cd(1);
161  h1Ch->DrawCopy();
162  c->cd(2);
163  hAllCh->DrawCopy();
164  c->cd(3);
165  h2D->DrawCopy("colz");
166 
167  c->cd(4);
168  hFFT1Ch->DrawCopy();
169  c->cd(5);
170  hFFTAllCh->DrawCopy();
171  c->cd(6);
172  hFFT2D->DrawCopy("colz");
173 
174  c->cd(7);
175  hPed1Ch->DrawCopy();
176  c->cd(8);
177  hPedAllCh->DrawCopy();
178  c->cd(9);
179  hPed2D->DrawCopy("colz");
180 
181  c2->cd(1);
182  h1Ev1Ch->DrawCopy();
183  c2->cd(2);
184  h1EvAllCh->DrawCopy();
185  c2->cd(3);
186  h1Ev2D->DrawCopy("colz");
187 
188  c2->cd(4);
189  h1EvFFT1Ch->DrawCopy();
190  c2->cd(5);
191  h1EvFFTAllCh->DrawCopy();
192  c2->cd(6);
193  h1EvFFT2D->DrawCopy("colz");
194 
195  c2->cd(7);
196  h1EvPed1Ch->DrawCopy();
197  c2->cd(8);
198  h1EvPedAllCh->DrawCopy();
199  c2->cd(9);
200  h1EvPed2D->DrawCopy("colz");
201 
202  TFile* file = new TFile("noise1Ev.root","recreate");
203 
204  h1Ev1Ch->Write();
205  h1EvAllCh->Write();
206  h1Ev2D->Write();
207 
208  h1EvFFT1Ch->Write();
209  h1EvFFTAllCh->Write();
210  h1EvFFT2D->Write();
211 
212  h1EvPed1Ch->Write();
213  h1EvPedAllCh->Write();
214  h1EvPed2D->Write();
215 
216  file->Close();
217 
218  h1Ev1Ch->Delete();
219  h1EvAllCh->Delete();
220  h1Ev2D->Delete();
221  h1EvFFT1Ch->Delete();
222  h1EvFFTAllCh->Delete();
223  h1EvFFT2D->Delete();
224  h1EvPed1Ch->Delete();
225  h1EvPedAllCh->Delete();
226  h1EvPed2D->Delete();
227 
228 
229 
230  } else {
231  std::cout << "No Data Map for plane " << std::endl;
232  }
233 
234  c->Update();
235  c2->Update();
236 
237 
238 
239  int iscan;
240  printf("Event %ld done. Waiting for input to continue: ",nEvents);
241  if(scanf("%d",&iscan) == 1)
242  std::cout << "Advance to next event" << std::endl;
243  else
244  std::cout << "Scan all events" << std::endl;
245 
246 }
247 
249 {
250  // m_hmap[0] = new HarpoMap(0);
251  // m_hmap[1] = new HarpoMap(1);
252 
253  h1Ch = new TH1F("h1Ch",";time [tb]",508,0,508);
254  hAllCh = new TH1F("hAllCh",";time [tb]",508,0,508);
255  h2D = new TH2F("h2D",";time [tb]; channel",508,0,508,304,0,304);
256 
257  hFFT1Ch = new TH1F("hFFT1Ch",";freq []",508,0,508);
258  hFFTAllCh = new TH1F("hFFTAllCh",";freq []",508,0,508);
259  hFFT2D = new TH2F("hFFT2D",";freq []; channel",508,0,508,304,0,304);
260 
261  hPed1Ch = new TH1F("hPed1Ch",";pedesteal [ADC]",500,0,500);
262  hPedAllCh = new TH1F("hPedAllCh",";pedesteal [ADC]",500,0,500);
263  hPed2D = new TH2F("hPed2D",";pedesteal [ADC]; channel",500,0,500,304,0,304);
264 
265 }
266 
268 {
269 
270 
271  TCanvas *c = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("cHarpo");
272  if(c){
273  c->SaveAs("noise.png");
274  c->SaveAs("noise.C");
275  }
276  // TString * hstFile = gHConfig->GetHistFile();
277  // if ( hstFile == NULL ) {
278  // std::cout << gHConfig->GetProgramName() << " " <<
279  // "No Hist File name gived, use default" << std::endl;
280  // hstFile = new TString("harpohst.root");
281  // }
282 
283  TFile* file = new TFile("noise.root","recreate");
284 
285  h1Ch->Write();
286  hAllCh->Write();
287  h2D->Write();
288 
289  hFFT1Ch->Write();
290  hFFTAllCh->Write();
291  hFFT2D->Write();
292 
293  hPed1Ch->Write();
294  hPedAllCh->Write();
295  hPed2D->Write();
296 
297  file->Close();
298 
299  // TFile *hf = new TFile(hstFile->Data(),"RECREATE");
300  // hf->Close();
301 }
302 
303 
304 
305 
306 
307 
308 
309 // Double_t HarpoFeminosAnalyseNoiseMonitor::TruncMean(TArrayD* vect, Double_t tl, Double_t th)
310 // {
311 // //
312 // // Calculates the truncated mean mean of the non zero value in vect
313 // // The truncation is done on the 100*tl% lowest and 100*(1-th) highest value
314 // //
315 // // debug(2,"truncated mean");
316 
317 
318 // Int_t size = vect->GetSize();
319 // Double_t truncMean = 0;
320 // Int_t* index = new Int_t[size];
321 // TMath::Sort(size,vect->GetArray(),index,kFALSE);
322 // Int_t t = 0, tLow, tHigh;
323 
324 
325 // while(vect->At(index[t])<10&&t<size-1) t++;
326 // tLow = TMath::FloorNint(t + (size - t)*tl);
327 // tHigh = TMath::FloorNint(t + (size - t)*th);
328 
329 // for(Int_t tind = tLow; tind<tHigh; tind++) truncMean += vect->At(index[tind]);
330 // if(tHigh-tLow) {
331 // return truncMean/(tHigh-tLow);
332 // }
333 // else return 0;
334 // }
335 
336 
337 
#define NALL
Definition: HarpoDccMap.h:15
Double_t GetData(Int_t i, Int_t j)
Set/Get Data Cell.
Definition: HarpoDccMap.cxx:84
virtual HarpoDccMap * GetMap()
unpacked dcc data The class contains the data map for DCC or Feminos The data is stored as a 2D TMatr...
Definition: HarpoDccMap.h:29
A virtual class which define intrafece between HARPO Reader and Event Analysis code.
#define NADC
Definition: HarpoDccMap.h:18
HarpoDetEvent * anaEvt