Splitting of the QA maker into simulation and reconstruction dependent parts (Yves)
[u/mrichter/AliRoot.git] / TRD / AliTRDQADataMakerRec.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //  Produces the data needed to calculate the quality assurance.          //
21 //  All data must be mergeable objects.                                   //
22 //                                                                        //
23 //  Author:                                                               //
24 //    Sylwester Radomski (radomski@physi.uni-heidelberg.de)               //
25 //                                                                        //
26 ////////////////////////////////////////////////////////////////////////////
27
28 // --- ROOT system ---
29 #include <TClonesArray.h>
30 #include <TFile.h> 
31 #include <TH1D.h> 
32 #include <TH2D.h>
33 #include <TH3D.h>
34 #include <TProfile.h>
35 #include <TF1.h>
36 #include <TCanvas.h>
37
38 // --- AliRoot header files ---
39 #include "AliESDEvent.h"
40 #include "AliLog.h"
41 #include "AliTRDcluster.h"
42 #include "AliTRDQADataMakerRec.h"
43 #include "AliTRDgeometry.h"
44 #include "AliTRDdataArrayI.h"
45 #include "AliTRDrawStreamTB.h"
46
47 #include "AliQAChecker.h"
48
49 ClassImp(AliTRDQADataMakerRec)
50
51 //____________________________________________________________________________ 
52   AliTRDQADataMakerRec::AliTRDQADataMakerRec() : 
53   AliQADataMakerRec(AliQA::GetDetName(AliQA::kTRD), "TRD Quality Assurance Data Maker")
54 {
55   //
56   // Default constructor
57 }
58
59 //____________________________________________________________________________ 
60 AliTRDQADataMakerRec::AliTRDQADataMakerRec(const AliTRDQADataMakerRec& qadm) :
61   AliQADataMakerRec()
62 {
63   //
64   // Copy constructor 
65   //
66
67   SetName((const char*)qadm.GetName()) ; 
68   SetTitle((const char*)qadm.GetTitle()); 
69
70 }
71
72 //__________________________________________________________________
73 AliTRDQADataMakerRec& AliTRDQADataMakerRec::operator=(const AliTRDQADataMakerRec& qadm)
74 {
75   //
76   // Equal operator.
77   //
78
79   this->~AliTRDQADataMakerRec();
80   new(this) AliTRDQADataMakerRec(qadm);
81   return *this;
82
83 }
84
85 //____________________________________________________________________________ 
86 void AliTRDQADataMakerRec::EndOfDetectorCycle(AliQA::TASKINDEX task, TObjArray * list)
87 {
88   //
89   // Detector specific actions at end of cycle
90   //
91
92   //AliInfo(Form("EndOfCycle", "Fitting RecPoints %d", task));
93
94  
95   if (task == AliQA::kRECPOINTS) {
96
97     //list->Print();
98     
99     // Rec points full chambers
100     if (((TH2D*)list->At(1))->GetEntries() > 1e4) {
101       for (Int_t i=0; i<540; i++) {
102         
103         TH1D *h = ((TH2D*)list->At(1))->ProjectionY(Form("qaTRD_recPoints_amp_%d",i), i+1, i+1);
104         if (h->GetSum() < 100) continue; // chamber not present
105         
106         h->Fit("landau", "q0", "goff", 10, 180);
107         TF1 *fit = h->GetFunction("landau");
108         ((TH1D*)list->At(12))->Fill(fit->GetParameter(1));
109         ((TH1D*)list->At(13))->Fill(fit->GetParameter(2));
110         delete h;      
111       }
112     }
113     
114     
115     if (((TH2D*)list->At(10))->GetEntries() > 1e5) {
116       for (Int_t i=0; i<540; i++) {
117         
118         TH1D *test = ((TH3D*)list->At(10))->ProjectionZ(Form("ampTime_%d",i), i+1, i+1, 0, 35);     
119         if (test->GetSum() < 100) continue;
120         
121         //AliInfo(Form("fitting det = %d", i));
122         
123         for(Int_t j=0; j<35; j++) {
124           
125           TH1D *h =  ((TH3D*)list->At(10))->ProjectionZ(Form("ampTime_%d",i), i+1, i+1, j+1, j+1);     
126           if (h->GetSum() < 50) continue;
127           
128           h->Fit("landau", "q0", "goff", 10, 180);
129           TF1 *fit = h->GetFunction("landau");
130           
131           Int_t sm = i/18;
132           Int_t det = i%18;
133           TH2D *h2 = (TH2D*)list->At(14+sm);
134           Int_t bin = h2->FindBin(det,j);
135           // printf("%d %d %d\n", det, j, bin);
136           h2->SetBinContent(bin, fit->GetParameter(1));
137         }
138       }
139     }
140   }
141   
142   // call the checker
143   AliQAChecker::Instance()->Run(AliQA::kTRD, task, list) ;    
144
145
146 }
147
148 //____________________________________________________________________________ 
149 void AliTRDQADataMakerRec::InitESDs()
150 {
151   //
152   // Create ESDs histograms in ESDs subdir
153   //
154
155   const Int_t kNhist = 19;
156   TH1 *hist[kNhist];
157   Int_t histoCounter = -1 ;
158
159   hist[++histoCounter] = new TH1D("qaTRD_esd_ntracks", ":Number of tracks", 300, -0.5, 299.5);
160   hist[++histoCounter] = new TH1D("qaTRD_esd_sector", ":Sector", 18, -0.5, 17.7);
161   hist[++histoCounter] = new TH1D("qaTRD_esd_bits", ";Bits", 64, -0.5, 63.5);
162
163   const Int_t knbits = 6;
164   const char *suf[knbits] = {"TPCi", "TPCo", "TPCz", "TRDo", "TRDr", "TRDz"};
165
166   for(Int_t i=0; i<knbits; i++) {
167     hist[++histoCounter] = new TH1D(Form("qaTRD_esd_pt%s",suf[i]), ";p_{T} (GeV/c);", 50, 0, 10);
168     hist[++histoCounter] = new TH1D(Form("qaTRD_esd_trdz%s", suf[i]), ";z (cm)", 200, -400, 400); 
169   }
170
171   hist[++histoCounter] = new TH1D("qaTRD_esd_clsTRDo", "TRDo;number of clusters", 130, -0.5, 129.5);;
172   hist[++histoCounter] = new TH1D("qaTRD_esd_clsTRDr", "TRDr;number of clusters", 130, -0.5, 129.5);;
173   hist[++histoCounter] = new TH1D("qaTRD_esd_clsTRDz", "TRDz;number of clusters", 130, -0.5, 129.5);;
174   //hist[++histoCounter] = new TH1D("qaTRD_esd_clsRatio", ";cluster ratio", 100, 0., 1.3);;
175
176   hist[++histoCounter] = new TH2D("qaTRD_esd_sigMom", ";momentum (GeV/c);signal", 100, 0, 5, 200, 0, 1e3);
177
178   for(Int_t i=0; i<=histoCounter; i++) {
179     //hist[i]->Sumw2();
180     Add2ESDsList(hist[i], i);
181   }
182
183 }
184
185 //____________________________________________________________________________ 
186 void AliTRDQADataMakerRec::InitRecPoints()
187 {
188   //
189   // Create Reconstructed Points histograms in RecPoints subdir
190   //
191
192   const Int_t kNhist = 14 + 18;
193   TH1 *hist[kNhist];
194
195   hist[0] = new TH1D("qaTRD_recPoints_det", ";Detector ID of the cluster", 540, -0.5, 539.5);
196   hist[1] = new TH2D("qaTRD_recPoints_amp", ";Amplitude", 540, -0.5, 539, 200, -0.5, 199.5);
197   hist[2] = new TH1D("qaTRD_recPoints_npad", ";Number of Pads", 12, -0.5, 11.5);
198
199   hist[3] = new TH1D("qaTRD_recPoints_dist2", ";residuals [2pad]", 100, -1, 1);
200   hist[4] = new TH1D("qaTRD_recPoints_dist3", ";residuals [3pad]", 100, -1, 1);
201   hist[5] = new TH1D("qaTRD_recPoints_dist4", ";residuals [4pad]", 100, -1, 1);
202   hist[6] = new TH1D("qaTRD_recPoints_dist5", ";residuals [5pad]", 100, -1, 1);
203
204   hist[7] = new TH2D("qaTRD_recPoints_rowCol", ";row;col", 16, -0.5, 15.5, 145, -0.5, 144.5);
205   hist[8] = new TH1D("qaTRD_recPoints_time", ";time bin", 35, -0.5, 34.5);
206   hist[9] = new TH1D("qaTRD_recPoints_nCls", ";number of clusters", 500, -0.5, 499.5);
207
208   hist[10] = new TH3D("qaTRD_recPoints_sigTime", ";chamber;time bin;signal", 
209                       540, -0.5, 539.5, 35, -0.5, 34.5, 100, 0, 200);
210   hist[11] = new TProfile("qaTRD_recPoints_prf", ";distance;center of gravity"
211                          , 120, -0.6, 0.6, -1.2, 1.2, "");
212
213   hist[12] = new TH1D("qaTRD_recPoints_ampMPV", ";amplitude MPV", 100, 0, 100);
214   hist[13] = new TH1D("qaTRD_recPoints_ampSigma", ";amplitude Sigma", 100, 0, 100); 
215
216   for(Int_t i=0; i<18; i++) {
217     hist[14+i] = new TH2D(Form("qaTRD_recPoints_sigTime_sm%d",i), Form("sm%d;det;time bin"), 
218                         30, -0.5, 29.5, 35, -0.5, 34.5);
219     hist[14+i]->SetMinimum(20);
220     hist[14+i]->SetMaximum(40);
221   }
222
223   for(Int_t i=0; i<kNhist; i++) {
224     //hist[i]->Sumw2();
225     Add2RecPointsList(hist[i], i);
226   }
227
228 }
229
230 //____________________________________________________________________________ 
231 void AliTRDQADataMakerRec::InitRaws()
232 {
233   //
234   // create Raws histograms in Raws subdir
235   //
236
237   const Int_t kSM = 18;
238   //const Int_t kNCh = 540;
239   const Int_t kNhist = 4+kSM;
240   TH1D *hist[kNhist];
241
242   // four histograms to be published
243   hist[0] = new TH1D("qaTRD_raws_det", ";detector", 540, -0.5, 539.5);
244   hist[1] = new TH1D("qaTRD_raws_sig", ";signal", 100, -0.5, 99.5);
245   hist[2] = new TH1D("qaTRD_raws_timeBin", ";time bin", 40, -0.5, 39.5); 
246   hist[3] = new TH1D("qaTRD_raws_smId", ";supermodule", 18, -0.5, 17.5);
247   //
248   
249   // one double per MCM (not published)
250   const Int_t kNMCM = 30 * 8 * 16;
251   for(Int_t i=0; i<kSM; i++)
252     hist[4+i] = new TH1D(Form("qaTRD_raws_sm%d",i),"",kNMCM, -0.5, kNMCM-0.5); 
253   
254   // register
255   for(Int_t i=0; i<kNhist; i++) {
256     //hist[i]->Sumw2();
257     Add2RawsList(hist[i], i);
258   }
259
260 }
261
262 //____________________________________________________________________________
263 void AliTRDQADataMakerRec::MakeESDs(AliESDEvent * esd)
264 {
265   //
266   // Make QA data from ESDs
267   //
268
269   Int_t nTracks = esd->GetNumberOfTracks();
270   GetESDsData(0)->Fill(nTracks);
271
272   // track loop
273   for (Int_t i=0; i<nTracks; i++) {
274
275     AliESDtrack *track = esd->GetTrack(i);
276     const AliExternalTrackParam *paramOut = track->GetOuterParam();
277     const AliExternalTrackParam *paramIn = track->GetInnerParam();
278
279     // long track ..
280     if (!paramIn) continue;
281     if (!paramOut) continue;
282
283     // not a kink
284     if (track->GetKinkIndex(0) > 0) continue; 
285
286     Double_t extZ = GetExtZ(paramIn);
287     if (TMath::Abs(extZ) > 320) continue; // acceptance cut
288
289     // .. in the acceptance
290     Int_t sector = GetSector(paramOut->GetAlpha());
291     GetESDsData(1)->Fill(sector);
292
293     UInt_t u = 1;
294     UInt_t status = track->GetStatus();
295     for(Int_t bit=0; bit<32; bit++) 
296       if (u<<bit & status) GetESDsData(2)->Fill(bit);
297
298     const Int_t knbits = 6; 
299     Int_t bit[6] = {0,0,0,0,0,0};    
300     bit[0] = status & AliESDtrack::kTPCin;
301     bit[1] = status & AliESDtrack::kTPCout;
302     bit[2] = (status & AliESDtrack::kTPCout) && !(status & AliESDtrack::kTRDout);
303     bit[3] = status & AliESDtrack::kTRDout;
304     bit[4] = status & AliESDtrack::kTRDrefit;
305     bit[5] = (status & AliESDtrack::kTRDout) && !(status & AliESDtrack::kTRDrefit);
306
307     // transverse momentum
308     //const Double_t *val = paramOut->GetParameter(); // parameters at the Outer plane
309     Double_t pt = paramOut->Pt(); //1./TMath::Abs(val[4]);
310
311     for(Int_t b=0; b<knbits; b++) {
312       if (bit[b]) {
313         GetESDsData(2*b+3)->Fill(pt); 
314         GetESDsData(2*b+4)->Fill(extZ);
315       }
316     }
317
318     // clusters
319     for(Int_t b=0; b<3; b++) 
320       if (bit[3+b]) GetESDsData(b+15)->Fill(track->GetTRDncls());
321
322     // refitted only
323     if (!bit[4]) continue;
324
325     //fQuality->Fill(track->GetTRDQuality());
326     //fBudget->Fill(track->GetTRDBudget());
327     //fSignal->Fill(track->GetTRDsignal());
328         
329     GetESDsData(18)->Fill(track->GetP(), track->GetTRDsignal());
330
331     /*
332     // PID only
333     if (status & AliESDtrack::kTRDpid) {
334
335       for(Int_t l=0; l<6; l++) fTime->Fill(track->GetTRDTimBin(l));
336
337       // fill pid histograms
338       Double_t trdr0 = 0; //, tpcr0 = 0;
339       Int_t trdBestPid = 5; //, tpcBestPid = 5;  // charged
340       const Double_t kminPidValue = 0.9;
341
342       //Double_t pp[5];
343       //track->GetTPCpid(pp); // ESD inconsequence
344
345       for(Int_t pid=0; pid<5; pid++) {
346         
347         trdr0 += track->GetTRDpid(pid);
348         //tpcr0 += pp[pid];
349         
350         fTrdPID[pid]->Fill(track->GetTRDpid(pid));
351         //fTpcPID[pid]->Fill(pp[pid]);
352         
353         if (track->GetTRDpid(pid) > kminPidValue) trdBestPid = pid;
354         //if (pp[pid] > kminPidValue) tpcBestPid = pid;
355       }
356
357       fTrdPID[5]->Fill(trdr0); // check unitarity
358       fTrdSigMomPID[trdBestPid]->Fill(track->GetP(), track->GetTRDsignal());
359
360       //fTpcPID[5]->Fill(tpcr0); // check unitarity
361       //fTpcSigMomPID[tpcBestPid]->Fill(track->GetP(), track->GetTPCsignal());
362     }
363     */
364
365   }
366
367 }
368
369 //______________________________________________________________________________
370 Int_t AliTRDQADataMakerRec::GetSector(const Double_t alpha) const 
371 {
372   //
373   // Gets the sector number 
374   //
375
376   Double_t size = TMath::DegToRad() * 20.; // shall use TRDgeo
377   Int_t sector = (Int_t)((alpha + TMath::Pi())/size);
378   return sector;
379
380 }
381
382 //______________________________________________________________________________
383 Double_t AliTRDQADataMakerRec::GetExtZ(const AliExternalTrackParam *in) const 
384 {
385   //
386   // Returns the Z position at the entry to TRD
387   // using parameters from the TPC in
388   //
389
390   const Double_t kX0 = 300;
391
392   Double_t x = in->GetX();
393   const Double_t *par = in->GetParameter();
394   Double_t theta = par[3];
395   Double_t z = in->GetZ();
396
397   Double_t zz = z + (kX0-x) * TMath::Tan(theta);
398   return zz;
399
400 }
401
402 //____________________________________________________________________________
403 void AliTRDQADataMakerRec::MakeRaws(AliRawReader* rawReader)
404 {
405   //
406   // Makes QA data from raw data
407   //
408
409   // 157
410   // T9 -- T10
411
412   //const Int_t kSM = 18;
413   //const Int_t kROC = 30;
414   const Int_t kROB = 8;
415   //const Int_t kLayer = 6;
416   //const Int_t kStack = 5;
417   const Int_t kMCM = 16;
418   //  const Int_t kADC = 22;
419
420   AliTRDrawStreamTB *raw = new AliTRDrawStreamTB(rawReader);
421
422   raw->SetRawVersion(3);
423   raw->Init();
424
425   while (raw->Next()) {
426
427     GetRawsData(0)->Fill(raw->GetDet());
428
429     // possibly needs changes with the new reader !!
430     Int_t *sig = raw->GetSignals();
431     for(Int_t i=0; i<3; i++) GetRawsData(1)->Fill(sig[i]);
432     // ---
433
434     GetRawsData(2)->Fill(raw->GetTimeBin());
435
436     // calculate the index;
437     Int_t sm = raw->GetSM();
438     Int_t roc = raw->GetROC();
439     Int_t rob = raw->GetROB();
440     Int_t mcm = raw->GetMCM();
441     //Int_t adc = raw->GetADC();
442
443     //Int_t index = roc * (kROB*kMCM*kADC) + rob * (kMCM*kADC) + mcm * kADC + adc;
444     Int_t  index = roc * (kROB*kMCM) + rob * kMCM + mcm;
445     GetRawsData(3)->Fill(sm);
446     GetRawsData(4+sm)->Fill(index);
447   }
448 }
449
450 //____________________________________________________________________________
451 void AliTRDQADataMakerRec::MakeRecPoints(TTree * clustersTree)
452 {
453   //  
454   // Makes data from RecPoints
455   // 
456
457   Int_t nsize = Int_t(clustersTree->GetTotBytes() / (sizeof(AliTRDcluster))); 
458   TObjArray *clusterArray = new TObjArray(nsize+1000); 
459
460   TBranch *branch = clustersTree->GetBranch("TRDcluster");
461   if (!branch) {
462     AliError("Can't get the branch !");
463     return;
464   }
465   branch->SetAddress(&clusterArray); 
466
467   // Loop through all entries in the tree
468   Int_t nEntries   = (Int_t) clustersTree->GetEntries();
469   Int_t nbytes     = 0;
470   AliTRDcluster *c = 0;
471   Int_t nDet[540];
472   for (Int_t i=0; i<540; i++) nDet[i] = 0;
473
474   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
475
476     // Import the tree
477     nbytes += clustersTree->GetEvent(iEntry);  
478
479     // Get the number of points in the detector
480     Int_t nCluster = clusterArray->GetEntriesFast();  
481
482     // Loop through all TRD digits
483     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
484       c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
485
486       Int_t iDet = c->GetDetector();
487       nDet[iDet]++;
488       GetRecPointsData(0)->Fill(iDet);
489       GetRecPointsData(1)->Fill(iDet, c->GetQ());
490       GetRecPointsData(2)->Fill(c->GetNPads());
491       if (c->GetNPads() < 6)
492         GetRecPointsData(1+c->GetNPads())->Fill(c->GetCenter());
493
494       //if (c->GetPadTime() < 5)
495       ((TH2D*)GetRecPointsData(7))->Fill(c->GetPadRow(), c->GetPadCol());
496       GetRecPointsData(8)->Fill(c->GetPadTime());
497
498       ((TH3D*)GetRecPointsData(10))->Fill(iDet, c->GetPadTime(), c->GetQ());
499
500       // PRF for 2pad
501       //if (c->GetNPads() == 2) {
502       Short_t *sig = c->GetSignals();
503       Double_t frac = -10;
504
505       if (sig[0] == 0 && sig[1] == 0 && sig[2] == 0 && sig[5] == 0 && sig[6] == 0) 
506         frac = 1. * sig[4] / (sig[3] + sig[4]);
507
508       if (sig[0] == 0 && sig[1] == 0 && sig[4] == 0 && sig[5] == 0 && sig[6] == 0)
509         frac = -1. * sig[2] / (sig[2] + sig[3]);
510
511       if (frac > -10)  ((TProfile*)GetRecPointsData(11))->Fill(c->GetCenter(), frac);
512         
513       //}
514     }
515   }
516
517   for(Int_t i=0; i<540; i++) 
518     if (nDet[i] > 0) GetRecPointsData(9)->Fill(nDet[i]);
519
520   delete clusterArray;
521
522 }
523
524 //____________________________________________________________________________ 
525 void AliTRDQADataMakerRec::StartOfDetectorCycle()
526 {
527   //
528   // Detector specific actions at start of cycle
529   //
530
531 }
532
533 //__________________________________________________________________________
534 Int_t AliTRDQADataMakerRec::CheckPointer(TObject *obj, const char *name) 
535 {
536   //
537   // Checks initialization of pointers
538   //
539
540   if (!obj) AliWarning(Form("null pointer: %s", name));
541   return !!obj;
542
543 }