1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TRD/AliTRDcluster.h>
3 #include <AliRunLoader.h>
10 #include <TObjArray.h>
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
22 // Number of noisy detectors
23 const Int_t numNoisyDets = 10;
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};
29 Bool_t isBadCluster(Int_t det, Int_t Q, Int_t PadCol, Int_t PadRow)
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;
38 // Exclude noisy detectors
39 for (Int_t indDet = 0; indDet < numNoisyDets; indDet++)
41 if (det == exDets[indDet])
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;
51 // Cluster seems to be clean
56 //=====================================
58 void AliTRDanalyseRecPoints (const char *filename="galice.root")
60 gROOT->SetStyle("Plain");
61 gStyle->SetPalette(1);
63 // open run loader and load gAlice, kinematics and header
64 AliRunLoader* runLoader = AliRunLoader::Open(filename);
67 printf("Error: readKine:\nGetting run loader from file \"%s\" failed", filename);
71 runLoader->LoadHeader();
72 runLoader->LoadRecPoints("TRD");
73 TObjArray *module = new TObjArray();
75 runLoader->CdGAFile();
78 // Flag indicating whether the signal comes from a noisy detector
79 Bool_t fromNoisyDet = 0;
81 // Counts the number of clusters per event
84 // Threshold for the amplitude
85 const Int_t kQThreshold = 30;
87 // Selection criterions for clusters
88 const Int_t kThresholdLower = 30;
89 const Int_t kThresholdUpper = 500;
91 // Number of events passing the criterion
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);
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);
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);
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);
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);
117 TH2D *mQvsPads = new TH2D("QvsPads", ";Number of pads;Q;#Clusters", 15, -0.5, 14.5, 100, 0, 1000);
119 // There are 540 detectors
120 TH2D *mQvsDetector = new TH2D("QvsDet", ";Number of detector;Q;#Clusters", 540, -0.5, 539.5, 100, 0, 1000);
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);
130 // Index + 2 corresponds to the number of activated pads
133 for (Int_t ind = 0; ind < 14; ind++)
135 sprintf(name,"%dPads",ind + 2);
136 mPads[ind] = new TH1D(name, ";Charge", 200, 0, 200);
139 // Amplitude distribution will be displayed event by event for the first 100 events
140 TH1D *mChrgEBE = new TH1D("chrgEBE", ";Charge", 200, 0, 200);
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);
145 // Keep track of the number of clusters per chamber with Q < kQThreshold
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);
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);
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);
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);
166 #ifdef EVENT_BY_EVENT
167 TCanvas* eventc = new TCanvas("eventc","Amplitude distribution for event 0",700,500);
171 #ifdef EVENT_BY_EVENT2
172 TCanvas* sTrc = new TCanvas("sTrc","Track - Supermodules",800,800);
179 int nEvents = runLoader->GetNumberOfEvents();
181 AliTRDcluster *cls = 0;
183 for(Int_t ev = 0; ev < nEvents - 1; ev++)
187 TTree *tree = runLoader->GetTreeR("TRD", 0);
188 tree->SetBranchAddress("TRDcluster", &module);
190 int N = tree->GetEntries();
192 // Check number of clusters
193 for(Int_t ind = 0; ind < N; ind++)
196 Int_t m = module->GetEntries();
198 for (Int_t j = 0; j < m; j++)
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++;
205 // Uncomment the following lines AND comment the next if-clause to display the histograms for the
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());
219 // Process event, if number of clusters is in the correct range
220 if (clusters >= kThresholdLower && clusters <= kThresholdUpper)
225 for (Int_t i = 0; i < N; i++)
228 int m = module->GetEntries();
234 for(Int_t j = 0; j < m; j++)
239 if (cls != 0) delete cls;
240 cls = (AliTRDcluster*)module->At(j);
242 // Filter the clusters of the events, too
243 if (cls->GetQ() < 20) continue;
245 mColTime->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
247 if (cls->GetQ() <= kQThreshold)
249 mClusCham->Fill(cls->GetDetector());
252 // Exclude noisy detectors
253 for (Int_t indDet = 0; indDet < numNoisyDets; indDet++)
255 if (cls->GetDetector() == exDets[indDet])
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;
269 mColRow->Fill(cls->GetPadCol(), cls->GetPadRow());
272 //if (cls->GetPadRow() == 15 || (cls->GetPadCol() >= 134 && cls->GetPadRow() <= 11) || (cls->GetPadCol() >= 119 && cls->GetPadCol() <= 123 && cls->GetPadRow() == 11)) continue;
274 if ((cls->GetPadCol() >= 134 && cls->GetPadRow() <= 5) || (cls->GetPadCol() == 141 && cls->GetPadRow() >= 10)) continue;
276 mColRowFiltered->Fill(cls->GetPadCol(), cls->GetPadRow());
279 mColTimeFiltered->Fill(cls->GetLocalTimeBin(), cls->GetPadCol());
281 if (cls->GetQ() <= kQThreshold)
283 mClusChamEx->Fill(cls->GetDetector());
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());
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());
301 mQvsPads->Fill(cls->GetNPads(),cls->GetQ());
302 mQvsDetector->Fill(cls->GetDetector(),cls->GetQ());
303 mQvsSModId->Fill((Int_t)(cls->GetDetector() / 30.),cls->GetQ());
305 Int_t superModule = (Int_t)(cls->GetDetector() / 30.);
307 // Position in supermodule = cls->GetDetector() - superModule * 30
308 // There are 6 layers per supermodule:
309 Int_t layer = (cls->GetDetector() - superModule * 30) % 6;
313 mSMod0->Fill(cls->GetQ());
314 mTrSMod0->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
317 mSMod8->Fill(cls->GetQ());
318 mTrSMod8->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
321 mSMod9->Fill(cls->GetQ());
322 mTrSMod9->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
325 mSMod17->Fill(cls->GetQ());
326 mTrSMod17->Fill(cls->GetPadCol(), cls->GetLocalTimeBin() + layer * 20);
332 // Only consider events with maximal 15 pads here
333 if (cls->GetNPads() - 2 < 14 && cls->GetNPads() - 2 >= 0)
335 mPads[cls->GetNPads() - 2]->Fill(cls->GetQ());
338 mChrgEBE->Fill(cls->GetQ());
343 #ifdef EVENT_BY_EVENT
344 // Event by event display for the first 100 events
347 sprintf(title,"Charge distribution for event %d", ev);
348 eventc->SetTitle(title);
356 #ifdef EVENT_BY_EVENT2
357 // Event by event display
358 sprintf(title2,"Supermodule tracks for event %d", ev);
359 sTrc->SetTitle(title2);
361 mTrSMod0->Draw("colz");
363 mTrSMod8->Draw("colz");
365 mTrSMod9->Draw("colz");
367 mTrSMod17->Draw("colz");
377 mClusters->Fill(clusters);
380 runLoader->GetNextEvent();
397 new TCanvas("chrgc", "Charge(no pad-filter)", 700, 500);
400 TF1* lFit = new TF1("lFit", "landau", 20, 200);
401 cout << "\n\nFit (no pad-filter):\n";
402 mChrg->Fit(lFit, "", "", 20, 200);
404 new TCanvas("chrg2c", "Charge(pads<9)", 700, 500);
407 TF1* lFit2 = new TF1("lFit2", "landau", 20, 200);
408 cout << "\n\nFit (pads<9):\n";
409 mChrg2->Fit(lFit2, "", "", 20, 200);
411 new TCanvas("chrg3c", "Charge(2<pads<9)", 700, 500);
414 TF1* lFit3 = new TF1("lFit3", "landau", 20, 200);
415 cout << "\n\nFit (2<pads<9):\n";
416 mChrg3->Fit(lFit3, "", "", 20, 200);
444 mPlaneY->Draw("colz");
449 mQvsPads->Draw("colz");
453 mQvsDetector->Draw("colz");
455 for (Int_t num = 30; num < 540; num += 30)
457 // Draw vertical lines every 30 bins to sketch the layer modules
458 tl = new TLine(num - 0.5, 0, num - 0.5, 1000);
464 mQvsSModId->Draw("colz");
466 TCanvas* sc = new TCanvas("sc","Supermodules", 800, 800);
479 TCanvas* padc = new TCanvas("padc","Pads", 1000, 1000);
481 for (Int_t ind = 0; ind < 14; ind++)
490 TCanvas* clusChamc = new TCanvas("clusChamc", "Clusters against detector with threshold", 800, 800);
491 clusChamc->Divide(1,2);
497 TCanvas* noisyc = new TCanvas("noisyc", "Col/Row of noisy detectors");
501 mColRow->Draw("colz");
504 mColRowFiltered->Draw("colz");
506 TCanvas* ColTimec = new TCanvas("ColTimec", "Col/Time");
507 ColTimec->Divide(1,2);
510 mColTime->Draw("colz");
513 mColTimeFiltered->Draw("colz");
515 cout << numPassed << " out of " << nEvents << " events passed the criterion (number of 'good' clusters >= ";
516 cout << kThresholdLower << " && clusters <= " << kThresholdUpper << ").\n";