next50 trigger mask in AliHLTGlobalEsdConverterComponent
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDanalyseRecPoints.C
CommitLineData
68252c5c 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
23const 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
27const Int_t exDets[numNoisyDets] = {25, 511, 515, 538, 523, 249, 294, 27, 24, 247};
28
29Bool_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
58void 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}