next50 trigger mask in AliHLTGlobalEsdConverterComponent
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDanalyseRecPoints.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TRD/AliTRDcluster.h>
3 #include <AliRunLoader.h>
4 #include <TCanvas.h>
5 #include <TF1.h>
6 #include <TH1.h>
7 #include <TH2.h>
8 #include <TLine.h>
9 #include <TMath.h>
10 #include <TObjArray.h>
11 #include <TPad.h>
12 #include <TROOT.h>
13 #include <TStyle.h>
14 #include <TSystem.h>
15 #include <TTree.h>
16 #endif
17
18 //#define EVENT_BY_EVENT      // Show amplitude event-by-event
19 //#define EVENT_BY_EVENT2     // Show tracks event-by-event
20 //#define PRINT_OTHER         // Display various additional histograms
21
22 // Number of noisy detectors
23 const Int_t numNoisyDets = 10;
24
25 // Holds the detector number of detectors that produce too much noise and
26 // are therefore excluded from the calculation
27 const Int_t exDets[numNoisyDets] = {25, 511, 515, 538, 523, 249, 294, 27, 24, 247};   
28
29 Bool_t isBadCluster(Int_t det, Int_t Q, Int_t PadCol, Int_t PadRow)
30 {  
31   if (Q < 20) return 1;
32
33   // Exclude bad cols
34   if (PadCol == 142 || PadCol == 136 || PadCol == 66 || PadCol == 65 || PadCol <= 2) return 1;
35   //if (PadCol == 142 || PadCol == 66 || PadCol == 65 || PadCol <= 2) return 1;
36   //if (PadCol == 142 || PadCol <= 2) return 1;
37   
38   // Exclude noisy detectors
39   for (Int_t indDet = 0; indDet < numNoisyDets; indDet++)
40   {
41     if (det == exDets[indDet])
42     {
43       // Exclude several regions of bad detectors
44       if ((PadCol >= 134 && PadRow <= 5) || (PadCol == 141 && PadRow >= 10)) return 1;
45       //if (PadRow == 15 || (PadCol >= 134 && PadRow <= 11) || (PadCol >= 119 && PadCol <= 123 && PadRow == 11)) return 1; 
46
47       break;
48     }
49   }
50   
51   // Cluster seems to be clean
52   return 0;
53 }
54
55
56 //=====================================
57
58 void AliTRDanalyseRecPoints (const char *filename="galice.root") 
59 {
60   gROOT->SetStyle("Plain");
61   gStyle->SetPalette(1);
62
63   // open run loader and load gAlice, kinematics and header
64   AliRunLoader* runLoader = AliRunLoader::Open(filename);
65   if (!runLoader) 
66   {
67     printf("Error: readKine:\nGetting run loader from file \"%s\" failed", filename);
68     return;
69   }
70
71   runLoader->LoadHeader();
72   runLoader->LoadRecPoints("TRD");
73   TObjArray *module = new TObjArray();
74
75   runLoader->CdGAFile();
76   
77
78   // Flag indicating whether the signal comes from a noisy detector
79   Bool_t fromNoisyDet = 0;
80
81   // Counts the number of clusters per event
82   Int_t clusters = 0;
83
84   // Threshold for the amplitude
85   const Int_t kQThreshold = 30;
86
87   // Selection criterions for clusters
88   const Int_t kThresholdLower = 30;
89   const Int_t kThresholdUpper = 500;
90
91   // Number of events passing the criterion
92   Int_t numPassed = 0;
93
94   // hist
95 #ifdef PRINT_OTHER
96   TH1D *mNCls = new TH1D("NCls", ";number of clusters", 400, 0, 400);
97   TH1D *mResY = new TH1D("resY", ";sigma Y (cm)", 100, 0, 1);
98   TH1D *mResZ = new TH1D("resZ", ";sigma Z (cm)", 100, 1, 3);
99
100   TH1D *mNPads = new TH1D("NPads", ";number of pads", 10, -0.5, 9.5);
101   TH1D *mCenter = new TH1D("center", ";center", 100, -1, 1);
102   TH1D *mTime = new TH1D("timeBin", ";time bin", 25, -0.5, 24.5);
103   TH1D *mPlane = new TH1D("plane", ";plane", 6, -0.5, 5.5);
104
105   TH1D *mX = new TH1D("x", "; X (cm)", 100, 0, 3.5);
106   TH1D *mY = new TH1D("y", "; y / chamber width (%)", 100, -1, 1);
107   TH1D *mZ = new TH1D("z", "; Z (cm)", 200, -400, 400);
108
109   TH2D *mPlaneY = new TH2D("planeY", ";plane;Y (cm);#Clusters", 6, -0.5, 5.5, 100, -60, 60);
110   TH2D *mYX = new TH2D("yx", ";Y (cm);X (cm);#Clusters", 100, -60, 60, 120, 280, 400);
111 #endif
112
113   TH1D *mChrg = new TH1D("chrg", ";Charge(no pad-filter)", 200, 0, 200);
114   TH1D *mChrg2 = new TH1D("chrg2", ";Charge(pads<9)", 200, 0, 200);
115   TH1D *mChrg3 = new TH1D("chrg3", ";Charge(2<pads<9)", 200, 0, 200);
116   
117   TH2D *mQvsPads = new TH2D("QvsPads", ";Number of pads;Q;#Clusters", 15, -0.5, 14.5, 100, 0, 1000);
118
119   // There are 540 detectors
120   TH2D *mQvsDetector = new TH2D("QvsDet", ";Number of detector;Q;#Clusters", 540, -0.5, 539.5, 100, 0, 1000);
121
122   // Supermoduleid=int(detectorNum / 30)
123   TH2D *mQvsSModId = new TH2D("QvsSModId", ";Supermodule id;Q;#Clusters", 18, -0.5, 17.5, 100, 0, 1000);
124   // Histograms for the already installed supermodules
125   TH1D *mSMod0 = new TH1D("SMod0", ";Charge", 200, 0, 200);
126   TH1D *mSMod8 = new TH1D("SMod8", ";Charge", 200, 0, 200);
127   TH1D *mSMod9 = new TH1D("SMod9", ";Charge", 200, 0, 200);
128   TH1D *mSMod17 = new TH1D("SMod17", ";Charge", 200, 0, 200);
129
130   // Index + 2 corresponds to the number of activated pads
131   TH1D *mPads[14];
132   char name[8];
133   for (Int_t ind = 0; ind < 14; ind++)
134   {
135     sprintf(name,"%dPads",ind + 2);
136     mPads[ind] = new TH1D(name, ";Charge", 200, 0, 200);
137   }
138
139   // Amplitude distribution will be displayed event by event for the first 100 events
140   TH1D *mChrgEBE = new TH1D("chrgEBE", ";Charge", 200, 0, 200);
141
142   // Keep track of the number of clusters per event
143   TH1D *mClusters = new TH1D("Clusters", ";Number of clusters per event", 300, -0.5, 299.5);
144
145   // Keep track of the number of clusters per chamber with Q < kQThreshold
146   char clusTitle[55];
147   sprintf(clusTitle,";Number of clusters with Q < %d vs detector", kQThreshold);
148   TH1D *mClusCham = new TH1D("ClusCham", clusTitle, 540, -0.5, 539.5);
149   // Same but this time noisy detectors excluded
150   sprintf(clusTitle,";Number of clusters with Q < %d vs detector", kQThreshold);
151   TH1D *mClusChamEx = new TH1D("ClusChamEx", clusTitle, 540, -0.5, 539.5);
152
153   TH2D *mColRow = new TH2D("ColRow", ";Column;Row;#Clusters", 144, -0.5, 143.5, 18, -0.5, 17.5);
154   TH2D *mColRowFiltered = new TH2D("ColRowFiltered", ";Column;Row;#Clusters", 144, -0.5, 143.5, 18, -0.5, 17.5);
155
156   TH2D *mColTime = new TH2D("ColTime", ";TimeBin;Column", 18, -0.5, 17.5, 144, -0.5, 143.5);
157   TH2D *mColTimeFiltered = new TH2D("ColTimeFiltered", ";TimeBin;Column;#Clusters", 18, -0.5, 17.5, 144, -0.5, 143.5);
158
159   // YX histograms for the already installed supermodules
160   TH2D *mTrSMod0 = new TH2D("TrSMod0", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
161   TH2D *mTrSMod8 = new TH2D("TrSMod8", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
162   TH2D *mTrSMod9 = new TH2D("TrSMod9", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
163   TH2D *mTrSMod17 = new TH2D("TrSMod17", ";Col;time+layer*20;#Clusters", 144, -0.5, 143.5, 120, -0.5, 119.5);
164
165
166 #ifdef EVENT_BY_EVENT  
167   TCanvas* eventc = new TCanvas("eventc","Amplitude distribution for event 0",700,500);
168   char title[40];
169 #endif
170
171 #ifdef EVENT_BY_EVENT2
172   TCanvas* sTrc = new TCanvas("sTrc","Track -  Supermodules",800,800);
173   sTrc->SetLogz();
174   sTrc->Divide(2,2);
175   char title2[40];
176 #endif
177
178  
179   int nEvents = runLoader->GetNumberOfEvents();
180   
181   AliTRDcluster *cls = 0;
182
183   for(Int_t ev = 0; ev < nEvents - 1; ev++) 
184   {
185     clusters = 0;    
186
187     TTree *tree = runLoader->GetTreeR("TRD", 0);
188     tree->SetBranchAddress("TRDcluster", &module);
189
190     int N = tree->GetEntries();
191
192     // Check number of clusters
193     for(Int_t ind = 0; ind < N; ind++)
194     {
195       tree->GetEntry(ind);
196       Int_t m = module->GetEntries();
197   
198       for (Int_t j = 0; j < m; j++)
199       {
200         if (cls != 0) delete cls;
201         cls = (AliTRDcluster*)module->At(j);
202         if (!isBadCluster(cls->GetDetector(), (Int_t)cls->GetQ(), cls->GetPadCol(), cls->GetPadRow())) clusters++;
203         else
204               {
205           // Uncomment the following lines AND comment the next if-clause to display the histograms for the 
206           // EXCLUDED clusters
207
208                 //mQvsPads->Fill(cls->GetNPads(),cls->GetQ());
209                 //mQvsDetector->Fill(cls->GetDetector(),cls->GetQ());
210                 //mQvsSModId->Fill((Int_t)(cls->GetDetector() / 30.),cls->GetQ());      
211                 //mColTime->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
212           //mColRow->Fill(cls->GetPadCol(), cls->GetPadRow());    
213                 //mClusCham->Fill(cls->GetDetector());
214           //mYX->Fill(cls->GetY(), cls->GetX());
215         }
216       }  
217     }
218
219     // Process event, if number of clusters is in the correct range
220     if (clusters >= kThresholdLower && clusters <= kThresholdUpper)
221     {
222       clusters = 0;
223       numPassed++;
224
225       for (Int_t i = 0; i < N; i++) 
226       {
227               tree->GetEntry(i);
228               int m = module->GetEntries();
229
230 #ifdef PRINT_OTHER
231               mNCls->Fill(m);
232 #endif
233
234               for(Int_t j = 0; j < m; j++) 
235               {
236           fromNoisyDet = 0;
237                 clusters++;
238
239           if (cls != 0) delete cls;
240                 cls = (AliTRDcluster*)module->At(j);
241               
242                 // Filter the clusters of the events, too
243           if (cls->GetQ() < 20) continue;
244            
245                 mColTime->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
246
247           if (cls->GetQ() <= kQThreshold) 
248           {
249             mClusCham->Fill(cls->GetDetector());
250                 }
251           
252           // Exclude noisy detectors
253           for (Int_t indDet = 0; indDet < numNoisyDets; indDet++)
254                 {
255                   if (cls->GetDetector() == exDets[indDet])
256                   {
257                   fromNoisyDet = 1;
258                     break;
259                   }
260                 }
261
262           // Exclude I
263           if (cls->GetPadCol() == 142 || cls->GetPadCol() == 136 || cls->GetPadCol() == 66 || cls->GetPadCol() == 65 || cls->GetPadCol() <= 2) continue;
264           //if (cls->GetPadCol() == 142 || cls->GetPadCol() == 66 || cls->GetPadCol() == 65 || cls->GetPadCol() <= 2) continue;
265           //if (cls->GetPadCol() == 142 || cls->GetPadCol() <= 2) continue;
266
267                 if (fromNoisyDet) 
268                 {
269                   mColRow->Fill(cls->GetPadCol(), cls->GetPadRow());    
270  
271                   // Exclude II
272                   //if (cls->GetPadRow() == 15 || (cls->GetPadCol() >= 134 && cls->GetPadRow() <= 11) || (cls->GetPadCol() >= 119 && cls->GetPadCol() <= 123 && cls->GetPadRow() == 11)) continue; 
273             
274             if ((cls->GetPadCol() >= 134 && cls->GetPadRow() <= 5) || (cls->GetPadCol() == 141 && cls->GetPadRow() >= 10)) continue; 
275             
276             mColRowFiltered->Fill(cls->GetPadCol(), cls->GetPadRow());
277                 }
278
279           mColTimeFiltered->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
280
281           if (cls->GetQ() <= kQThreshold) 
282           {
283             mClusChamEx->Fill(cls->GetDetector());
284                 }
285                 
286 #ifdef PRINT_OTHER
287                 mX->Fill(cls->GetX());
288                 mZ->Fill(cls->GetZ());
289           mYX->Fill(cls->GetY(), cls->GetX());
290                 mResY->Fill(TMath::Sqrt(cls->GetSigmaY2()));
291                 mResZ->Fill(TMath::Sqrt(cls->GetSigmaZ2()));
292           mNPads->Fill(cls->GetNPads());
293           mCenter->Fill(cls->GetCenter());
294                 mTime->Fill(cls->GetLocalTimeBin());
295 #endif
296
297                 if (cls->GetNPads() > 2 && cls->GetNPads() < 9) mChrg3->Fill(cls->GetQ());
298           if (cls->GetNPads() < 9) mChrg2->Fill(cls->GetQ());
299           mChrg->Fill(cls->GetQ());              
300                 
301                 mQvsPads->Fill(cls->GetNPads(),cls->GetQ());
302                 mQvsDetector->Fill(cls->GetDetector(),cls->GetQ());
303                 mQvsSModId->Fill((Int_t)(cls->GetDetector() / 30.),cls->GetQ());        
304
305           Int_t superModule = (Int_t)(cls->GetDetector() / 30.);
306
307           // Position in supermodule = cls->GetDetector() - superModule * 30
308           // There are 6 layers per supermodule:
309           Int_t layer = (cls->GetDetector() - superModule * 30) % 6;
310                 switch (superModule)
311                 {
312                   case 0:
313                     mSMod0->Fill(cls->GetQ());
314                     mTrSMod0->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
315                     break;
316                   case 8:
317                     mSMod8->Fill(cls->GetQ());
318                     mTrSMod8->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
319                     break;
320                   case 9:
321                     mSMod9->Fill(cls->GetQ());
322                     mTrSMod9->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
323                     break;
324                   case 17:
325                     mSMod17->Fill(cls->GetQ());
326                     mTrSMod17->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
327                     break;
328                   default:
329                     break;
330                 }
331
332                 // Only consider events with maximal 15 pads here
333                 if (cls->GetNPads() - 2 < 14 && cls->GetNPads() - 2 >= 0)
334                 {
335                   mPads[cls->GetNPads() - 2]->Fill(cls->GetQ());
336                 }
337
338                 mChrgEBE->Fill(cls->GetQ());
339               }
340       }
341       
342
343 #ifdef EVENT_BY_EVENT  
344       // Event by event display for the first 100 events
345       if (numPassed < 100)
346       {
347         sprintf(title,"Charge distribution for event %d", ev);
348         eventc->SetTitle(title);
349         mChrgEBE->Draw();
350         eventc->Update();
351         gSystem->Sleep(400);
352         mChrgEBE->Reset();
353       }
354 #endif
355
356 #ifdef EVENT_BY_EVENT2
357       // Event by event display 
358       sprintf(title2,"Supermodule tracks for event %d", ev);
359       sTrc->SetTitle(title2);
360       sTrc->cd(1);
361       mTrSMod0->Draw("colz");
362       sTrc->cd(2);
363       mTrSMod8->Draw("colz");
364       sTrc->cd(3);
365       mTrSMod9->Draw("colz");
366       sTrc->cd(4);
367       mTrSMod17->Draw("colz");
368       sTrc->Update();
369       gSystem->Sleep(400);
370       mTrSMod0->Reset();
371       mTrSMod8->Reset();
372       mTrSMod9->Reset();
373       mTrSMod17->Reset();
374 #endif
375
376
377       mClusters->Fill(clusters);    
378     }
379      
380     runLoader->GetNextEvent();
381   }
382
383 #ifdef PRINT_OTHER
384   new TCanvas();
385   gPad->SetLogy();
386   mNCls->Draw();
387
388   new TCanvas();
389   gPad->SetLogy();
390   mResY->Draw();
391
392   new TCanvas();
393   gPad->SetLogy();
394   mResZ->Draw();
395 #endif
396
397   new TCanvas("chrgc", "Charge(no pad-filter)", 700, 500);
398   mChrg->Draw();
399   // Try landau-fit
400   TF1* lFit = new TF1("lFit", "landau", 20, 200);
401   cout << "\n\nFit (no pad-filter):\n";
402   mChrg->Fit(lFit, "", "", 20, 200);
403
404   new TCanvas("chrg2c", "Charge(pads<9)", 700, 500);
405   mChrg2->Draw();
406   // Try landau-fit
407   TF1* lFit2 = new TF1("lFit2", "landau", 20, 200);
408   cout << "\n\nFit (pads<9):\n";
409   mChrg2->Fit(lFit2, "", "", 20, 200);
410
411   new TCanvas("chrg3c", "Charge(2<pads<9)", 700, 500);
412   mChrg3->Draw();
413   // Try landau-fit
414   TF1* lFit3 = new TF1("lFit3", "landau", 20, 200);
415   cout << "\n\nFit (2<pads<9):\n";
416   mChrg3->Fit(lFit3, "", "", 20, 200);
417
418 #ifdef PRINT_OTHER
419   new TCanvas();
420   mNPads->Draw();
421
422   new TCanvas();
423   mTime->Draw();
424
425   new TCanvas();
426   mX->Draw();
427
428   new TCanvas();
429   mY->Draw();
430
431   new TCanvas();
432   mZ->Draw();
433
434   new TCanvas();
435   mCenter->Draw();
436
437   new TCanvas();
438   mYX->Draw("colz");
439
440   new TCanvas();
441   mPlane->Draw();
442
443   new TCanvas();
444   mPlaneY->Draw("colz");
445 #endif
446
447   new TCanvas();
448   gPad->SetLogz();
449   mQvsPads->Draw("colz");  
450
451   new TCanvas();
452   gPad->SetLogz();
453   mQvsDetector->Draw("colz");
454   TLine* tl = 0;
455   for (Int_t num = 30; num < 540; num += 30)
456   {
457     // Draw vertical lines every 30 bins to sketch the layer modules
458     tl = new TLine(num - 0.5, 0, num - 0.5, 1000);
459     tl->Draw();
460   }
461
462   new TCanvas();
463   gPad->SetLogz();
464   mQvsSModId->Draw("colz");
465
466   TCanvas* sc = new TCanvas("sc","Supermodules", 800, 800);
467   sc->SetLogz();
468   sc->Divide(2,2);
469   sc->cd(1);
470   mSMod0->Draw();
471   sc->cd(2);
472   mSMod8->Draw();
473   sc->cd(3);
474   mSMod9->Draw();
475   sc->cd(4);
476   mSMod17->Draw();
477  
478
479   TCanvas* padc = new TCanvas("padc","Pads", 1000, 1000);
480   padc->Divide(4,4);
481   for (Int_t ind = 0; ind < 14; ind++)
482   {
483     padc->cd(ind + 1);
484     mPads[ind]->Draw();
485   }
486
487   new TCanvas();
488   mClusters->Draw();
489
490   TCanvas* clusChamc = new TCanvas("clusChamc", "Clusters against detector with threshold", 800, 800);
491   clusChamc->Divide(1,2);
492   clusChamc->cd(1);
493   mClusCham->Draw();
494   clusChamc->cd(2);
495   mClusChamEx->Draw();
496
497   TCanvas* noisyc = new TCanvas("noisyc", "Col/Row of noisy detectors");
498   noisyc->Divide(1,2);
499   noisyc->cd(1);
500   gPad->SetLogz();
501   mColRow->Draw("colz");
502   noisyc->cd(2);
503   gPad->SetLogz();
504   mColRowFiltered->Draw("colz");
505
506   TCanvas* ColTimec = new TCanvas("ColTimec", "Col/Time");
507   ColTimec->Divide(1,2);
508   ColTimec->cd(1);
509   gPad->SetLogz();
510   mColTime->Draw("colz");
511   ColTimec->cd(2);
512   gPad->SetLogz();
513   mColTimeFiltered->Draw("colz");
514
515   cout << numPassed << " out of " << nEvents << " events passed the criterion (number of 'good' clusters >= ";
516   cout << kThresholdLower << " && clusters <= " << kThresholdUpper << ").\n";
517 }