HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoAnalyseNoiseMonitor.cxx
Go to the documentation of this file.
1 //
2 // File HarpoAnalyseNoiseMonitor.cxx
3 //
12 #include "HarpoAnalyseNoiseMonitor.h"
13 #include "HarpoConfig.h"
14 #include "HarpoDetSet.h"
15 //#include "HarpoDetSet.h"
16 #include "HarpoEvent.h"
17 
18 #include "TROOT.h"
19 #include "TStyle.h"
20 #include "TApplication.h"
21 #include "TFile.h"
22 #include "TCanvas.h"
23 #include "TLatex.h"
24 #include "TGraph.h"
25 #include "TH2F.h"
26 #include "TF1.h"
27 #include "TLinearFitter.h"
28 #include "TMath.h"
29 
30 #include <cstdlib>
31 #include <cstring>
32 #include <cassert>
33 #include <iostream>
34 
35 ClassImp(HarpoAnalyseNoiseMonitor)
36 
37 void HarpoAnalyseNoiseMonitor::print()
38 {
39  UInt_t nd; // number of detectors
40  HarpoEventHeader *hdr = fEvt->GetHeader();
41 
42  assert(hdr != NULL);
43  hdr->print();
44 
45  for (nd = 0; nd < gkNDetectors; nd++) {
46  // if (fEvt->isdataExist(nd)) {
47  HarpoDccMap *plane = fEvt->GetDccMap(nd);
48  if (plane != NULL )plane->print();
49  }
50 }
51 
52 void HarpoAnalyseNoiseMonitor::process()
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","Traces",0,0,1000,800);
64  c->Clear();
65  c->Divide(2,4);//,0.001,0.001);
66  }
67 
68  TCanvas *c2 = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("cHarpo2");
69  if(!c2){
70  c2 = new TCanvas("cHarpo2","Parameters",0,0,1000,800);
71  c2->Clear();
72  c2->Divide(2,5);//,0.001,0.001);
73  }
74 
75  gStyle->SetOptStat(0);
76 
77  TH1F* hTev[2];
78  //hTev[0] = (TH1F *)gROOT->GetListOfCanvases()->FindObject("hTev_X");
79  //if(!hTev[0])
80  hTev[0] = new TH1F("hTev_X","",511,0,511);
81  //hTev[1] = (TH1F *)gROOT->GetListOfCanvases()->FindObject("hTev_Y");
82  //if(!hTev[1])
83  hTev[1] = new TH1F("hTev_Y","",511,0,511);
84  // hTev[0]->Reset();
85  // hTev[1]->Reset();
86 
87  TH1F* hQtotTev[2];
88  Int_t nbins = 1000;
89  Int_t Qmax = 10000;
90  //hQtotTev[0] = (TH1F *)gROOT->GetListOfCanvases()->FindObject("hQtotTev_X");
91  //if(!hQtotTev[0])
92  hQtotTev[0] = new TH1F("hQtotTev_X","",nbins,0,Qmax);
93  //hQtotTev[1] = (TH1F *)gROOT->GetListOfCanvases()->FindObject("hQtotTev_Y");
94  //if(!hQtotTev[1])
95  hQtotTev[1] = new TH1F("hQtotTev_Y","",nbins,0,Qmax);
96 
97  hTev[0]->SetLineColor(kBlue);
98  hTev[1]->SetLineColor(kRed);
99  hQtotTev[0]->SetLineColor(kBlue);
100  hQtotTev[1]->SetLineColor(kRed);
101 
102  //TFile* fileFFT = new TFile("testFFT/test.root","recreate");
103 
104 
105  // Double_t qtot[2] = {-1,-1};
106  // Double_t qtb[NADC][2];
107  //Int_t npix = 0, nch = 0, ntb = 0;
108  for(UInt_t ndet=0;ndet<gkNDetectors;ndet++) {
109  if (gHDetSet->isExist(ndet)) {
110  HarpoDccMap * m = fEvt->GetDccMap(ndet);
111  std::cout << "Event " << nEvents << ", detector " << ndet << " " << m << std::endl;
112  if ( m != NULL ) {
113  // for(j=0;j<NCHAN;j++) // Channels
114  Double_t qTot = 0;
115  TH2F* hFull = new TH2F("hFull",";time [25ns];channel [1mm]",511,0,511,304,0,304);
116  TH2F* hMinusTB = new TH2F("hMinusTB",";time [25ns];channel [1mm]",511,0,511,304,0,304);
117  TH2F* hMinusCh = new TH2F("hMinusCh",";time [25ns];channel [1mm]",511,0,511,304,0,304);
118  TH2F* hMinusAll = new TH2F("hMinusAll",";time [25ns];channel [1mm]",511,0,511,304,0,304);
119  //TH2F* hMinusAll = new TH2F("hMinusAll",";time [25ns];channel [1mm]",304,0,304,1000,-200,200);
120  TH1F* hPar0 = new TH1F("hPar0",";par0",1000,-100,100);
121  TH1F* hPar1 = new TH1F("hPar1",";par1",1000,-1,1);
122  TH1F* hPar2 = new TH1F("hPar2",";par2",1000,-0.01,0.01);
123  TH1F* hQ = new TH1F("hQ",";rms",1000,0,50);
124  TH1F* hQcut = new TH1F("hQcut",";rms",1000,0,50);
125  TH1F* hNoiseTB = new TH1F("hNoiseTB",";noiseTB",1000,0,1000);
126  // c->cd(7+ndet);
127  // hMinusAll->DrawCopy("colz");
128  // TH1F* hNoise[304];
129  // Double_t noiseTB[NADC];
130  // Double_t noiseParCh[NADC][3];
131  // Double_t noiseCh[NALL];
132  for(Int_t i=0;i<NADC;i++){ //Time bins
133  Double_t qTotT = 0;
134  TArrayD* valuesTB = new TArrayD(NALL);
135  for(Int_t j=0;j<NALL;j++){ // Channels
136  // if(i==0 && j<288) hNoise[j] = new TH1F(Form("hNoise_%d_%d",ndet,j),";time [25ns];q",511,0,511);
137  Double_t q = m->GetData(j,i);
138  if(ndet==1)
139  q = m->GetData(2*Int_t(j/2)+1-j%2,i);
140  hFull->SetBinContent(i+1,j+1,q);
141  valuesTB->AddAt(q,j);
142  qTotT += q;
143  qTot += q;
144  }
145  //noiseTB[i] = TruncMean(valuesTB,0.16,0.84);
146  Double_t noiseTB = TruncMean(valuesTB,0.16,0.84);
147  hNoiseTB->Fill(noiseTB);
148  TGraph* gMinusCh = new TGraph();
149  for(Int_t j=0;j<NALL;j++){ // Channels
150  // if(i==0 && j<288) hNoise[j] = new TH1F(Form("hNoise_%d_%d",ndet,j),";time [25ns];q",511,0,511);
151  Double_t q = m->GetData(j,i);
152  if(ndet==1)
153  q = m->GetData(2*Int_t(j/2)+1-j%2,i);
154  q -= noiseTB;
155  hMinusTB->SetBinContent(i+1,j+1,q);
156  gMinusCh->SetPoint(gMinusCh->GetN(),j,q);
157  }
158  // c->cd(7+ndet);
159  // gMinusCh->SetMarkerStyle(2);
160  // gMinusCh->SetMarkerColor(51+Int_t(i*49./(NADC-1)));
161  // gMinusCh->Draw("Psame");
162 
163  TLinearFitter* fFitCh = new TLinearFitter(1,"pol2","D");
164  fFitCh->AssignData(gMinusCh->GetN(), 1, gMinusCh->GetX(), gMinusCh->GetY());
165  fFitCh->Eval();
166  //fFitCh->EvalRobust();
167  // TF1* fFitCh = new TF1("fFitCh","pol2",0,288);
168  // gMinusCh->Fit(fFitCh,"q0");
169  Double_t mean = 0, rms = 0;
170  Int_t nnoise = 0;
171  for(Int_t j=0;j<NALL;j++){ // Channels
172  Double_t q = m->GetData(j,i);
173  if(ndet==1)
174  q = m->GetData(2*Int_t(j/2)+1-j%2,i);
175  q -= noiseTB;
176  q -= fFitCh->GetParameter(0) + j*fFitCh->GetParameter(1) + j*j*fFitCh->GetParameter(2);
177  hMinusCh->SetBinContent(i+1,j+1,q);
178  if(q<100){
179  mean += q;
180  rms += q*q;
181  nnoise++;
182  }
183  }
184  rms = TMath::Sqrt(rms/nnoise - mean*mean/nnoise/nnoise);
185  for(Int_t j=0;j<NALL;j++){ // Channels
186  Double_t q = hMinusCh->GetBinContent(i+1,j+1);
187  hQ->Fill(q);
188  if(q>3*rms){
189  hQcut->Fill(q);
190  hMinusAll->SetBinContent(i+1,j+1,q);
191  }else
192  hMinusAll->SetBinContent(i+1,j+1,0);
193  }
194 
195  //Info("process","Fit param %f %f %f",fFitCh->GetParameter(0),fFitCh->GetParameter(1),fFitCh->GetParameter(2));
196 
197  hPar0->Fill(fFitCh->GetParameter(0));
198  hPar1->Fill(fFitCh->GetParameter(1));
199  hPar2->Fill(fFitCh->GetParameter(2));
200 
201  // for(Int_t k=0;k<3;k++)
202  // noiseParCh[i][k] = fFitCh->GetParameter(k);
203  fFitCh->Delete();
204 
205  }
206  // for(Int_t j=0;j<NALL;j++){ // Channels
207  // TArrayD* valuesCh = new TArrayD(NADC);
208  // for(Int_t i=0;i<NADC;i++){ //Time bins
209  // Double_t q = m->GetData(j,i);
210  // if(ndet==1)
211  // q = m->GetData(2*Int_t(j/2)+1-j%2,i);
212  // q -= noiseTB[i];
213  // valuesCh->AddAt(q,i);
214  // hMinusTB->SetBinContent(i+1,j+1,q);
215  // }
216  // noiseCh[j] = TruncMean(valuesCh,0.16,0.84);
217  // //Info("process","noiseCh %d: %f",j,noiseCh[j]);
218  // for(Int_t i=0;i<NADC;i++){ //Time bins
219  // Double_t q = m->GetData(j,i);
220  // if(ndet==1)
221  // q = m->GetData(2*Int_t(j/2)+1-j%2,i);
222  // q -= noiseTB[i];
223  // q -= noiseParCh[i][0] + j*noiseParCh[i][1] + j*j*noiseParCh[i][2];
224  // // q -= noiseCh[j];
225  // hMinusCh->SetBinContent(i+1,j+1,q);
226  // }
227  // }
228 
229  if(qTot>0)
230  hQtot[ndet]->Fill(qTot);
231  // qtot[ndet] = qTot;
232 
233  c->cd(1+ndet);
234  hFull->DrawCopy("colz");
235  c->cd(3+ndet);
236  hMinusTB->DrawCopy("colz");
237  c->cd(5+ndet);
238  hMinusCh->DrawCopy("colz");
239  c->cd(7+ndet);
240  hMinusAll->DrawCopy("colz");
241 
242 
243  c2->cd(1+ndet);
244  hNoiseTB->DrawCopy();
245  c2->cd(3+ndet);
246  hPar0->DrawCopy();
247  c2->cd(5+ndet);
248  hPar1->DrawCopy();
249  c2->cd(7+ndet);
250  hPar2->DrawCopy();
251  c2->cd(9+ndet);
252  hQ->DrawCopy();
253  hQcut->SetLineColor(kRed);
254  hQcut->DrawCopy("same");
255 
256  hFull->Delete();
257  hMinusTB->Delete();
258  hMinusCh->Delete();
259  hMinusAll->Delete();
260 
261  hNoiseTB->Delete();
262  hPar0->Delete();
263  hPar1->Delete();
264  hPar2->Delete();
265  hQ->Delete();
266  hQcut->Delete();
267  } else {
268  std::cout << "No Data Map for plane "<< ndet << std::endl;
269  }
270  }
271  }
272 
273  // c->SaveAs(Form("plots/checks/noise_%d.C",nEvents));
274  // c->SaveAs(Form("plots/checks/noisePar_%d.C",nEvents));
275 
276 
277  //fileFFT->Close();
278 
279 
280  // c2->SaveAs("testY.C");
281 
282 
283  // for(i=0;i<NADC;i++)
284  // hQtotXY->Fill(qtb[i][0],qtb[i][1]);
285 
286  // if(npix>600)
287 
288  // c->cd(3);
289  // hTev[0]->DrawCopy();
290  // hTev[1]->DrawCopy("same");
291  // hTev[0]->Delete();
292  // hTev[1]->Delete();
293  // c->cd(4);
294  // hT[0]->Draw();
295  // hT[1]->Draw("same");
296 
297  // c->SaveAs("test.png");
298 
299 
300  // TLatex* latex = new TLatex();
301  // latex->SetTextFont(132);
302  // latex->SetTextAlign(12);
303  // latex->SetTextSize(0.07);
304  // c->cd(1);
305  // latex->DrawLatex(20,280,"X direction");
306  // c->cd(2);
307  // latex->DrawLatex(20,280,"Y direction");
308  // c->cd(0);
309  // latex->DrawLatex(0.4,0.95,Form("Run %d, event %d",gHConfig->GetRunNo(),nEvents));
310 
311 
312  // c2->cd(1);
313  // hQtotTev[0]->DrawCopy();
314  // hQtotTev[1]->DrawCopy("same");
315  // hQtotTev[0]->Delete();
316  // hQtotTev[1]->Delete();
317  // hQ[0]->Draw();
318  // hQ[1]->Draw("same");
319  // c2->GetPad(1)->SetLogy();
320  // c2->cd(2);
321  // hQtotT[0]->Draw();
322  // hQtotT[1]->Draw("same");
323  // c2->GetPad(2)->SetLogx();
324 
325  // c2->cd(3);
326  // hCh[0]->Draw();
327  // hCh[1]->Draw("same");
328  // c2->cd(4);
329  // c2->GetPad(4)->SetLogy();
330  // hQtot[0]->Draw();
331  // hQtot[1]->Draw("same");
332  // c2->GetPad(4)->SetLogx();
333 
334 
335  // // fEvt->print();
336  // Long_t ndet;
337 
338  // HarpoDccMap *harpodccmap;
339 
340  // for(ndet=0;ndet<2;ndet++)
341  // {
342  // if(gHDetSet->isExist(ndet)==0) continue;
343  // //
344  // harpodccmap = fEvt->GetDccMap(ndet);
345  // if(!harpodccmap) continue;
346  // //
347  // m_hmap[ndet]->FillRawMap(harpodccmap);
348  // m_hmap[ndet]->SubtractPedestals();
349  // }
350 
351  // char commentname[1000];
352  // sprintf(commentname,"Evt %d",fEvt->GetHeader()->GetEvtNo());
353 
354  // m_hmap[0]->ShowMaps(commentname);
355  // m_hmap[1]->ShowMaps(commentname);
356 
357  c->Update();
358  c2->Update();
359 
360 
361  int iscan;
362  printf("Event %ld done. Waiting for input to continue: ",nEvents);
363  if(scanf("%d",&iscan) == 1)
364  std::cout << "Advance to next event" << std::endl;
365  else
366  std::cout << "Scan all events" << std::endl;
367 
368 }
369 
370 void HarpoAnalyseNoiseMonitor::Init()
371 {
372  // m_hmap[0] = new HarpoMap(0);
373  // m_hmap[1] = new HarpoMap(1);
374 
375  hT[0] = new TH1F("hT_X",";time [tb]",511,0,511);
376  hT[1] = new TH1F("hT_Y",";time [tb]",511,0,511);
377 
378  hCh[0] = new TH1F("hCh_X",";channel",304,0,304);
379  hCh[1] = new TH1F("hCh_Y",";channel",304,0,304);
380 
381  hQ[0] = new TH1F("hQ_X",";Q_{pixel} [ADC]",4096,0,4096);
382  hQ[1] = new TH1F("hQ_Y",";Q_{pixel} [ADC]",4096,0,4096);
383 
384 
385  const Int_t nbins = 500;
386  const Int_t nbins2 = 500;
387  Double_t xmin = 100;
388  Double_t xmin2 = 1e4;
389  Double_t xmax = 2e5;
390  Double_t xmax2 = 1e8;
391  Double_t logxmin = TMath::Log(xmin);
392  Double_t logxmin2 = TMath::Log(xmin2);
393  Double_t logxmax = TMath::Log(xmax);
394  Double_t logxmax2 = TMath::Log(xmax2);
395  Double_t binwidth = (logxmax-logxmin)/nbins;
396  Double_t binwidth2 = (logxmax2-logxmin2)/nbins2;
397  Double_t xbins[nbins+1];
398  Double_t xbins2[nbins+1];
399  xbins[0] = xmin;
400  xbins2[0] = xmin2;
401  for (Int_t i=1;i<=nbins;i++)
402  xbins[i] = xmin + TMath::Exp(logxmin+i*binwidth);
403  for (Int_t i=1;i<=nbins2;i++)
404  xbins2[i] = xmin2 + TMath::Exp(logxmin2+i*binwidth2);
405 
406  hQtotT[0] = new TH1F("hQtotT_X",";Q_{cluster} [ADC]",nbins,xbins);
407  hQtotT[1] = new TH1F("hQtotT_Y",";Q_{cluster} [ADC]",nbins,xbins);
408 
409 
410  hQtot[0] = new TH1F("hQtot_X","Q_{event} [ADC]",nbins2,xbins2);
411  hQtot[1] = new TH1F("hQtot_Y","Q_{event} [ADC]",nbins2,xbins2);
412 
413  hCh[0]->SetLineColor(kBlue);
414  hT[0]->SetLineColor(kBlue);
415  hQtotT[0]->SetLineColor(kBlue);
416  hQtot[0]->SetLineColor(kBlue);
417  hQ[0]->SetLineColor(kBlue);
418  hCh[1]->SetLineColor(kRed);
419  hT[1]->SetLineColor(kRed);
420  hQtotT[1]->SetLineColor(kRed);
421  hQtot[1]->SetLineColor(kRed);
422  hQ[1]->SetLineColor(kRed);
423 
424 
425 }
426 
427 void HarpoAnalyseNoiseMonitor::Save(char * /* mode */)
428 {
429  // TString * hstFile = gHConfig->GetHistFile();
430  // if ( hstFile == NULL ) {
431  // std::cout << gHConfig->GetProgramName() << " " <<
432  // "No Hist File name gived, use default" << std::endl;
433  // hstFile = new TString("harpohst.root");
434  // }
435 
436  // TFile *hf = new TFile(hstFile->Data(),"RECREATE");
437  // hf->Close();
438 }
439 
440 
441 
442 
443 
444 
445 
446 Double_t HarpoAnalyseNoiseMonitor::TruncMean(TArrayD* vect, Double_t tl, Double_t th)
447 {
448  //
449  // Calculates the truncated mean mean of the non zero value in vect
450  // The truncation is done on the 100*tl% lowest and 100*(1-th) highest value
451  //
452  // debug(2,"truncated mean");
453 
454 
455  Int_t size = vect->GetSize();
456  Double_t truncMean = 0;
457  Int_t* index = new Int_t[size];
458  TMath::Sort(size,vect->GetArray(),index,kFALSE);
459  Int_t t = 0, tLow, tHigh;
460 
461 
462  while(vect->At(index[t])<10&&t<size-1) t++;
463  tLow = TMath::FloorNint(t + (size - t)*tl);
464  tHigh = TMath::FloorNint(t + (size - t)*th);
465 
466  for(Int_t tind = tLow; tind<tHigh; tind++) truncMean += vect->At(index[tind]);
467  if(tHigh-tLow) {
468  return truncMean/(tHigh-tLow);
469  }
470  else return 0;
471 }
472 
473 
474 
NALL
#define NALL
Definition: HarpoDccMap.h:15
HarpoAnalyseNoiseMonitor
Monitor noise in non-zerosuppressed data.
Definition: HarpoAnalyseNoiseMonitor.h:17
HarpoDccMap::GetData
Double_t GetData(Int_t i, Int_t j)
Set/Get Data Cell.
Definition: HarpoDccMap.cxx:84
HarpoDetSet::isExist
Bool_t isExist(ULong_t det)
Detecror date exist //! Number of Real Detectors.
Definition: HarpoDetSet.h:33
HarpoDccMap::print
virtual void print()
Definition: HarpoDccMap.cxx:124
HarpoEventHeader
FullEvent Header not scecific to the detectors The class is ....
Definition: HarpoEventHeader.h:27
HarpoEventHeader::print
virtual void print()
Definition: HarpoEventHeader.cxx:21
HarpoDccMap
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
NADC
#define NADC
Definition: HarpoDccMap.h:18
HarpoAnalyse::fEvt
HarpoEvent * fEvt
Definition: HarpoAnalyse.h:70
HarpoAnalyse::nEvents
ULong_t nEvents
Definition: HarpoAnalyse.h:75
HarpoAnalyseNoiseMonitor::TruncMean
Double_t TruncMean(TArrayD *vect, Double_t tl, Double_t th)
Redefine empty default.
Definition: HarpoAnalyseNoiseMonitor.cxx:446
gkNDetectors
const ULong_t gkNDetectors
Definition: HarpoDet.h:14
HarpoEvent::GetDccMap
HarpoDccMap * GetDccMap(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:108
gHDetSet
R__EXTERN HarpoDetSet * gHDetSet
Definition: HarpoDetSet.h:71
Generated by   doxygen 1.8.6