Adding Domenico Colella as responsible for SPD part in TRI pp
[u/mrichter/AliRoot.git] / TRD / TRDrec / 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 <TH1F.h> 
32 #include <TH2F.h>
33 #include <TH3F.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 "AliRawReader.h"
42 #include "AliTRDcluster.h"
43 #include "AliTRDQADataMakerRec.h"
44 #include "AliTRDgeometry.h"
45 #include "AliTRDrawStream.h"
46
47 #include "AliTRDdigitsManager.h"
48 #include "AliTRDSignalIndex.h"
49 #include "AliTRDarrayADC.h"
50
51 #include "AliQAChecker.h"
52
53 ClassImp(AliTRDQADataMakerRec)
54
55 //____________________________________________________________________________ 
56   AliTRDQADataMakerRec::AliTRDQADataMakerRec() : 
57   AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kTRD), "TRD Quality Assurance Data Maker")
58 {
59   //
60   // Default constructor
61 }
62
63 //____________________________________________________________________________ 
64 AliTRDQADataMakerRec::AliTRDQADataMakerRec(const AliTRDQADataMakerRec& qadm) :
65   AliQADataMakerRec()
66 {
67   //
68   // Copy constructor 
69   //
70
71   SetName((const char*)qadm.GetName()) ; 
72   SetTitle((const char*)qadm.GetTitle()); 
73
74 }
75
76 //__________________________________________________________________
77 AliTRDQADataMakerRec& AliTRDQADataMakerRec::operator=(const AliTRDQADataMakerRec& qadm)
78 {
79   //
80   // Equal operator.
81   //
82
83   this->~AliTRDQADataMakerRec();
84   new(this) AliTRDQADataMakerRec(qadm);
85   return *this;
86
87 }
88
89 //____________________________________________________________________________ 
90 void AliTRDQADataMakerRec::EndOfDetectorCycle(AliQAv1::TASKINDEX_t task, TObjArray ** list)
91 {
92   //
93   // Detector specific actions at end of cycle
94   //
95   //TStopwatch watch;
96   //watch.Start();
97   ResetEventTrigClasses();
98   AliDebug(AliQAv1::GetQADebugLevel(), "End of TRD cycle");
99   //
100   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {  // loop over event types
101     //
102     if (!IsValidEventSpecie(specie, list)) continue;
103     SetEventSpecie(AliRecoParam::ConvertIndex(specie));
104     //
105     for (int itc=-1;itc<GetNTrigClasses();itc++) { // RS: loop over eventual clones per trigger class
106       //
107       if (task == AliQAv1::kRAWS) {             
108         TH2 *mnCls   = (TH2*)GetRawsData(1, itc);
109         TH2 *mClsDet = (TH2*)GetRawsData(2, itc); 
110         if (mnCls && mClsDet) mClsDet->Divide(mnCls);
111       }
112       //      
113       if (task == AliQAv1::kRECPOINTS) {
114         //
115         TH1F * hist = new TH1F("fitHist", "", 200, -0.5, 199.5);
116         TObjArray& arrRP = *GetRecPointsDataOfTrigClass(itc); // RS Histos matching to trigger class
117         // fill detector map;
118         TH1* h0     = (TH1*) arrRP[0];
119         TH1* detMap = (TH1*) arrRP[87];
120         if (h0 && detMap) {
121           for(Int_t i = 0 ; i < 540 ; i++) {
122             Double_t v = h0->GetBinContent(i+1);
123             Int_t sm = i/30;
124             Int_t det = i%30;
125             Int_t bin = detMap->FindBin(sm, det);
126             detMap->SetBinContent(bin, v);
127           }
128         }
129         //        
130         // Rec points full chambers
131         TH2* h2tmp = (TH2*) arrRP[1];
132         TH1* h12   = (TH1*) arrRP[12];
133         TH1* h13   = (TH1*) arrRP[13];
134         if (h2tmp && h12 && h13) {
135           for (Int_t i = 0 ; i < 540 ; i++) {
136             //AliDebug(AliQAv1::GetQADebugLevel(), Form("I = %d", i));
137             hist->Reset();          
138             // project TH2F into TH1F 
139             for(Int_t b = 1 ; b < hist->GetXaxis()->GetNbins()-1 ; b++) {
140               Double_t xvalue = hist->GetBinCenter(b);
141               Int_t bin = h2tmp->FindBin(i,xvalue);
142               Double_t value =  h2tmp->GetBinContent(bin);
143               hist->SetBinContent(b, value);
144             }
145             //AliDebug(AliQAv1::GetQADebugLevel(), Form("Sum = %d %f\n", i, hist->GetSum()));
146             if (hist->GetSum() < 100) continue; // not enougth data in a chamber
147             //
148             hist->Fit("landau", "q0", "goff", 10, 180);
149             TF1 *fit = hist->GetFunction("landau");
150             h12->Fill(fit->GetParameter(1));
151             h13->Fill(fit->GetParameter(2));
152           }
153         }
154         //        
155         // time-bin by time-bin sm by sm
156         TH3* h3tmp = (TH3*) arrRP[10]; 
157         if (h3tmp) {
158           for(Int_t i=0; i<18; i++) { // loop over super-modules      
159             for(Int_t j=0; j<kTimeBin; j++) { // loop over time bins      
160               hist->Reset();
161               for(Int_t b = 1 ; b < hist->GetXaxis()->GetNbins()-1 ; b++) {
162                 Double_t xvalue = hist->GetBinCenter(b);
163                 Double_t svalue = 0.0;
164                 for(Int_t det = i*30 ; det < (i+1)*30 ; det++) { // loop over detectors
165                   Int_t bin = h3tmp->FindBin(det,j,xvalue);
166                   Double_t value =  h3tmp->GetBinContent(bin);
167                   svalue += value;
168                 }
169                 //AliDebug(AliQAv1::GetQADebugLevel(), Form("v = %f\n", value));
170                 hist->SetBinContent(b, svalue);
171               }
172               //
173               if (hist->GetSum() < 100) continue;
174               //
175               hist->Fit("landau", "q0", "goff", 10, 180);
176               TF1 *fit = hist->GetFunction("landau");
177               //
178               TH1 *hi = (TH1*)arrRP[14+18+i];
179               if (!hi) continue;
180               hi->SetMarkerStyle(20);
181               Int_t bin = hi->FindBin(j);
182               // printf("%d %d %d\n", det, j, bin);
183               //
184               Double_t value = TMath::Abs(fit->GetParameter(1));
185               Double_t error = TMath::Abs(fit->GetParError(1));
186               //
187               if (value/error < 3) continue; // insuficient statistics
188               //
189               hi->SetBinContent(bin, value);
190               hi->SetBinError(bin, error);      
191             } // j
192           } // i
193           //
194           // for numerical convergence
195           TF1 *form = new TF1("formLandau", "landau", 0, 200);
196           //
197           // time-bin by time-bin chamber by chamber
198           for (Int_t i=0; i<540; i++) {
199             for(Int_t j=0; j<kTimeBin; j++) {
200               //            
201               hist->Reset();
202               for(Int_t b = 1 ; b < hist->GetXaxis()->GetNbins()-1 ; b++) {
203                 Double_t xvalue = hist->GetBinCenter(b);
204                 Int_t bin = h3tmp->FindBin(i,j,xvalue);
205                 Double_t value =  h3tmp->GetBinContent(bin);
206                 //AliDebug(AliQAv1::GetQADebugLevel(), Form("v = %f\n", value));
207                 hist->SetBinContent(b, value);
208               }
209               //
210               if (hist->GetSum() < 100) continue;
211               
212               form->SetParameters(1000, 60, 20);
213               hist->Fit(form, "q0", "goff", 20, 180);
214               
215               Int_t sm = i/30;
216               Int_t det = i%30;
217               TH1 *hi = (TH1*)arrRP[14+sm];
218               if (!hi) continue;
219               Int_t bin = hi->FindBin(det,j);
220               // printf("%d %d %d\n", det, j, bin);
221               //
222               Double_t value = TMath::Abs(form->GetParameter(1));
223               Double_t error = TMath::Abs(form->GetParError(1));
224               //
225               if (value/error < 3) continue;
226               //
227               hi->SetBinContent(bin, value);
228               hi->SetBinError(bin, error);
229             } // j
230           } // i
231         } // h3tmp      
232         if (hist) delete hist;
233       } // RECPOINTS
234       //
235       //////////////////////////
236       // const Int_t knbits = 6;
237       // const char *suf[knbits] = {"TPCi", "TPCo", "TPCz", "TRDo", "TRDr", "TRDz"};
238       //const char *sufRatio[4] = {"TRDrTRDo", "TRDoTPCo", "TRDrTPCo", "TRDzTPCo"};
239       
240       if (task == AliQAv1::kESDS) {
241         //
242         const Int_t knRatio = 4;
243         const Int_t kN[knRatio] = {4,3,4,5};
244         const Int_t kD[knRatio] = {3,1,1,3}; 
245         //
246         TObjArray& arrES = *GetESDsDataOfTrigClass(itc); // RS Histos matching to trigger class
247         // create ratios
248         for(Int_t type = 0 ; type < 2 ; type++) {
249           for(Int_t i = 0 ; i < knRatio ; i++) {            
250             TH1 *ratio = (TH1*)arrES[19 + 2*i + type];
251             TH1 *histN = (TH1*)arrES[3 + 2*kN[i] + type];
252             TH1 *histD = (TH1*)arrES[3 + 2*kD[i] + type];
253             if ( ! (ratio && histN && histD) ) continue;
254             BuildRatio(ratio, histN, histD);
255             //ratio->Reset();
256             //ratio->Add(histN);
257             //ratio->Divide(histD);
258           }
259         }
260         // ratio for the fraction of electrons per stack
261         TH1 *histN = (TH1*)arrES[33];
262         TH1 *histD = (TH1*)arrES[32];
263         TH1 *ratio = (TH1*)arrES[34];
264         BuildRatio(ratio, histN, histD);
265       } // ESDs
266     } // RS: loop over eventual clones per trigger class
267   } // loop over species
268   // call the checker
269   AliQAChecker::Instance()->Run(AliQAv1::kTRD, task, list) ;    
270   //
271 }
272
273 //____________________________________________________________________________ 
274 void AliTRDQADataMakerRec::InitESDs()
275 {
276   //
277   // Create ESDs histograms in ESDs subdir
278   //
279   const Bool_t expert   = kTRUE ; 
280   const Bool_t image    = kTRUE ; 
281   
282   const Int_t kNhist = 36+5+4;
283
284   TH1 *hist[kNhist];
285   Int_t histoCounter = -1 ;
286
287   hist[++histoCounter] = new TH1F("qaTRD_esd_ntracks", "TRD esd ntracks;Number of tracks;Counts", 300, -0.5, 299.5);
288   hist[++histoCounter] = new TH1F("qaTRD_esd_sector", "TRD esd sector;Sector;Counts", 18, -0.5, 17.7);
289   hist[++histoCounter] = new TH1F("qaTRD_esd_bits", "TRD esd bits;Bits;Counts", 64, -0.5, 63.5);
290
291   const Int_t knbits = 6;
292   const char *suf[knbits] = {"TPCi", "TPCo", "TPCz", "TRDo", "TRDr", "TRDz"};
293
294   // histo = 3
295   for(Int_t i=0; i<knbits; i++) { 
296     hist[++histoCounter] = new TH1F(Form("qaTRD_esd_pt%s",suf[i]), Form("qaTRD_esd_pt%s;p_{T} (GeV/c);Counts",suf[i]), 100, 0, 10);
297     hist[++histoCounter] = new TH1F(Form("qaTRD_esd_trdz%s", suf[i]), ";z (cm)", 200, -400, 400); 
298   }
299
300   hist[++histoCounter] = new TH1F("qaTRD_esd_clsTRDo", "TRDo;number of clusters;Counts", 180, -0.5, 179.5);;
301   hist[++histoCounter] = new TH1F("qaTRD_esd_clsTRDr", "TRDr;number of clusters;Counts", 180, -0.5, 179.5);;
302   hist[++histoCounter] = new TH1F("qaTRD_esd_clsTRDz", "TRDz;number of clusters;Counts", 180, -0.5, 179.5);;
303   //hist[++histoCounter] = new TH1F("qaTRD_esd_clsRatio", ";cluster ratio", 100, 0., 1.3);;
304
305   hist[++histoCounter] = new TH2F("qaTRD_esd_sigMom", "TRD esd sig Mom;momentum (GeV/c);signal", 100, 0, 5, 200, 0, 1e3);
306
307   // end of cycle plots (hist 19)
308   const char *sufRatio[4] = {"TRDrTRDo", "TRDoTPCo", "TRDrTPCo", "TRDzTPCo"};
309
310   for(int i=0; i<4; i++) {
311     hist[++histoCounter] = new TH1F(Form("qaTRD_esd_pt%s",sufRatio[i]), 
312                                     Form("Efficiency in Pt %s;p_{T};eff", sufRatio[i]),
313                                     100, 0, 10);
314
315     hist[++histoCounter] = new TH1F(Form("qaTRD_esd_trdz%s",sufRatio[i]), 
316                                     Form("Efficiency in Z %s;z (cm);eff", sufRatio[i]),
317                                     200, -400, 400);
318   }
319
320   // 27 - 31
321   hist[27] = new TH1F("qaTRD_esd_quality", "TRD esd quality;quality;Counts", 120, 0, 12);
322   hist[28] = new TH1F("qaTRD_esd_budget", "TRD esd budget;NN;Counts", 110, -1000, 100);
323   hist[29] = new TH1F("qaTRD_esd_chi2", "TRD esd chi2;chi2;Counts", 200, 0, 100);
324   hist[30] = new TH1F("qaTRD_esd_timeBin", "TRD esd timeBin;time bin;Counts", 7, -0.5, 6.5);
325   hist[31] = new TH1F("qaTRD_esd_pidQuality", "pid Quality;quality;;Counts", 7, -0.5, 6.5);
326
327   // stack by stack electron identyfication
328   hist[32] = new TH1F("qaTRD_esd_tracksStack", "number of all tracks;stack;Counts", 90, -0.5, 89.5);
329   hist[33] = new TH1F("qaTRD_esd_electronStack", "number of electron tracks;stack;Counts", 90, -0.5, 89.5);
330   hist[34] = new TH1F("qaTRD_esd_elRatioStack", "fraction of electron tracks;stack;Counts", 90, -0.5, 89.5);
331   hist[35] = new TH1F("qaTRD_esd_thetaOut", ";tan(theta);", 100, -1, 1);
332   
333   const char *partType[5] = {"Electron", "Muon", "Pion", "Kaon", "Proton"}; 
334
335   for(Int_t i=0; i<AliPID::kSPECIES; i++)
336     hist[36+i] = new TH1F(Form("qaTRD_esd_pid%d",i),
337                           Form("%s;probability;Counts",partType[i]), 100, 0, 1);
338  
339   // dE/dX vs momentum in three regions
340   const char *zoneName[4] = {"total charge", "ampilification range", "plateau", "TR range"};
341  
342   // prepare the scale from 0.1 to 10 GeV
343   const Int_t nscalex= 50;
344   Double_t scalex[nscalex+1];
345   Double_t dd = (TMath::Log(10) - TMath::Log(0.5)) / nscalex;
346   for(Int_t ix=0; ix<nscalex+1; ix++) {
347     scalex[ix] = 0.5 * TMath::Exp(dd * ix);
348   }
349
350   const Int_t nscaley = 50;
351   Double_t scaley[nscaley+1];
352   for(Int_t iy=0; iy<nscaley+1; iy++) {
353     scaley[iy] = iy * (3e3/nscaley);
354   }
355     
356
357   for(Int_t i=0; i<4; i++) {
358     hist[41+i] = new TH2F(Form("qaTRD_esd_signalPzone_%d",i), 
359                           Form("%s;momentum (GeV/c);signal (a.u.)", zoneName[i]),
360                           nscalex, scalex, nscaley, scaley);
361   }
362
363   for(Int_t i=0; i<kNhist; i++) {
364     //hist[i]->Sumw2();
365     Add2ESDsList(hist[i], i, !expert, image);
366   }
367   //
368   ClonePerTrigClass(AliQAv1::kESDS); // this should be the last line
369 }
370
371 //____________________________________________________________________________ 
372 void AliTRDQADataMakerRec::InitRecPoints()
373 {
374   //
375   // Create Reconstructed Points histograms in RecPoints subdir
376   //
377   const Bool_t expert   = kTRUE ; 
378   const Bool_t image    = kTRUE ; 
379   
380   //printf("Helo from Init rec points\n");
381
382   const Int_t kNhist = 14 + 4 * 18 + 2 + 9;// + 540;
383   TH1 *hist[kNhist];
384
385   hist[0] = new TH1F("qaTRD_recPoints_det", "RRD recPoints det;Detector ID of the cluster;Counts", 540, -0.5, 539.5);
386   hist[1] = new TH2F("qaTRD_recPoints_amp", "TRD recPoints amp;Amplitude;??", 540, -0.5, 539, 200, -0.5, 199.5);
387   hist[2] = new TH1F("qaTRD_recPoints_npad", "TRD recPoints npad;Number of Pads;Counts", 12, -0.5, 11.5);
388
389   hist[3] = new TH1F("qaTRD_recPoints_dist2", "TRD recPoints dist2;residuals [2pad];Counts", 100, -1, 1);
390   hist[4] = new TH1F("qaTRD_recPoints_dist3", "TRD recPoints dist3;residuals [3pad];Counts", 100, -1, 1);
391   hist[5] = new TH1F("qaTRD_recPoints_dist4", "TRD recPoints dist4;residuals [4pad];Counts", 100, -1, 1);
392   hist[6] = new TH1F("qaTRD_recPoints_dist5", "TRD recPoints dist5;residuals [5pad];Counts", 100, -1, 1);
393
394   hist[7] = new TH2F("qaTRD_recPoints_rowCol", "TRDrecPointsrowCol;row;col", 16, -0.5, 15.5, 145, -0.5, 144.5);
395   hist[8] = new TH1F("qaTRD_recPoints_time", "TRDrecPoints time;time bin;Counts", kTimeBin, -0.5, kTimeBin-0.5);
396   hist[9] = new TH1F("qaTRD_recPoints_nCls", "TRD recPoints nCls;number of clusters;Counts", 500, -0.5, 499.5);
397
398   hist[10] = new TH3F("qaTRD_recPoints_sigTime", "TRD recPoints sigTime;chamber;time bin;signal", 
399                       540, -0.5, 539.5, kTimeBin, -0.5, kTimeBin-0.5, 200, -0.5, 199.5);
400   hist[11] = new TProfile("qaTRD_recPoints_prf", "TRD recPoints prf;distance;center of gravity;Counts"
401                          , 120, -0.6, 0.6, -1.2, 1.2, "");
402
403   hist[12] = new TH1F("qaTRD_recPoints_ampMPV", "TRD recPoints ampMPV;amplitude MPV;Counts", 150, 0, 150);
404   hist[13] = new TH1F("qaTRD_recPoints_ampSigma", "TRD recPoints ampSigma;amplitude Sigma;Counts", 200, 0, 200); 
405   
406   // chamber by chamber
407   for(Int_t i=0; i<18; i++) {
408     hist[14+i] = new TH2F(Form("qaTRD_recPoints_sigTime_sm%d",i), Form("sm%d;det;time bin",i), 
409                         30, -0.5, 29.5, kTimeBin, -0.5, kTimeBin-0.5);
410     hist[14+i]->SetMinimum(0);
411     hist[14+i]->SetMaximum(150);
412   }
413  
414   // time bin by time bin sm-by-sm
415   for(Int_t i=0; i<18; i++) {
416     hist[14+18+i] = new TH1F(Form("qaTRD_recPoints_sigTimeShape_sm%d", i), 
417                              Form("sm%d;time bin;signal",i),
418                              kTimeBin, -0.5, kTimeBin-0.5);
419
420     hist[14+18+i]->SetMaximum(150);    
421   }
422
423   // str = 50
424   for(Int_t i=0; i<18; i++) {
425     hist[50+i] = new TH1F(Form("qaTRD_recPoints_nCls_sm%d",i),
426                           Form("sm%d;time bin;number of clusters",i),
427                           kTimeBin, -0.5, kTimeBin-0.5);
428   }
429
430   // str = 68
431   for(Int_t i=0; i<18; i++) {
432     hist[68+i] = new TH1F(Form("qaTRD_recPoints_totalCharge_sm%d", i),
433                           Form("sm%d;time bin;total charge", i),
434                           kTimeBin, -0.5, kTimeBin-0.5);
435   }
436
437   hist[86] = new TH1F("qaTRD_recPoints_signal", "TRD recPoints signal;amplitude;Counts", 400, -0.5, 399.5);
438   hist[87] = new TH2F("qaTRD_recPoints_detMap", "TRD recPoints detMap;sm;chamber;Counts", 18, -0.5, 17.5, 30, -0.5, 29.5);
439
440   
441   // amplitude as a function of the pad size
442   for(Int_t i=0; i<9; i++) {
443     hist[88+i] = new TH1F(Form("qaTRD_recPoints_signalNpad_%d", i+2), Form("qaTRD_recPoints_signalNpad_%d;amplitude, ADC", i+2), 400, -0.5, 399.5); 
444   }
445   
446   // one 2D histogram per chamber
447   //  for(Int_t i=0; i<540; i++) {
448   //  hist[88+i] = new TH2F(Form("qaTRD_recPoints_map%d", i), ";col;row", 16, -0.5, 15.5, 144, -0.5, 143.5);
449   //}
450
451
452   for(Int_t i=0; i<kNhist; i++) {
453     //hist[i]->Sumw2();
454     Add2RecPointsList(hist[i], i, !expert, image);
455   }
456   //
457   ClonePerTrigClass(AliQAv1::kRECPOINTS); // this should be the last line
458 }
459
460 //____________________________________________________________________________ 
461 void AliTRDQADataMakerRec::InitRaws()
462 {
463   //
464   // create Raws histograms in Raws subdir
465   //
466   const Bool_t expert   = kTRUE ; 
467   const Bool_t saveCorr = kTRUE ; 
468   const Bool_t image    = kTRUE ; 
469   
470   AliInfo("Initialization of QA for Raw Data");
471   
472   const Int_t kNhist = 7;
473   TH1 *hist[kNhist];
474
475   hist[0] = new TH2F("qaTRD_raws_nADC","number of ADC channels;sector;detector", 18, -0.5, 17.5, 30, -0.5, 29.5);
476   hist[1] = new TH2F("qaTRD_raws_nCls", "number of clusters;sector;detector", 18, -0.5, 17.5, 30, -0.5, 29.5);
477   hist[2] = new TH2F("qaTRD_raws_meanSig", "mean signal;sector;detector", 18, -0.5, 17.5, 30, -0.5, 29.5);
478   
479   hist[3] = new TH1F("qaTRD_raws_ADC", "ADC amplitude;ADC counts", 100, -0.5, 99.5);
480   hist[4] = new TH1F("qaTRD_raws_Cls", "Cluster amplitude; ADC counts", 100, -0.5, 199.5);
481   hist[5] = new TH2F("qaTRD_raws_ClsTb", "Clusters vs Time Bin;time bin;amoplitude", 30, -0.5, 29.5, 200, -0.5, 199.5);
482   
483   hist[6] = new TH2F("qaTRD_raws_ClsAmpDet", ";detector;amplitude", 540, -0.5, 539.5, 100, 0, 200);
484   
485
486   /*
487     hist[0] = new TH2F("qaTRD_raws_DataVolume", ";Sector;Data Volume, kB", 18, -0.5, 17.5, 100, 0, 30);
488     hist[1] = new TH2F("qaTRD_raws_HC", "Data Headers;Sector;HC", 18, -0.5, 17.5, 60, -0.5, 59.5);
489     hist[2] = new TH2F("qaTRD_raws_LME", "Link Monitor Error;Sector;HC", 18, -0.5, 17.5, 60, -0.5, 59.5);
490     
491     hist[3] = new TH1F("qaTRD_rawd_cls", "Clusters amplitude;ADC counts", 100, 0, 200);
492     hist[4] = new TH2F("qaTRD_raws_clsTB", "amplitude - time bins;time bin;amplitude"
493     30, -0.5, 29.5, 100, 0, 200);
494     hist[5] = new TH2F("qaTRD_raws_clsSec", "amplitude in sectors;Sector;amplitude, ADCs"
495     18, -0.5, 17.5);
496   */
497
498
499   // register
500   for(Int_t i=0; i<kNhist; i++) {
501     //hist[i]->Sumw2();
502     Add2RawsList(hist[i], i, !expert, image, !saveCorr);
503   }
504   //
505   ClonePerTrigClass(AliQAv1::kRAWS); // this should be the last line
506 }
507
508 //____________________________________________________________________________
509 void AliTRDQADataMakerRec::MakeESDs(AliESDEvent * esd)
510 {
511   //
512   // Make QA data from ESDs
513   //
514   
515   Int_t nTracks = esd->GetNumberOfTracks();
516   FillESDsData(0,nTracks);
517
518   // track loop
519   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++) {
520
521     AliESDtrack *track = esd->GetTrack(iTrack);
522     const AliExternalTrackParam *paramOut = track->GetOuterParam();
523     const AliExternalTrackParam *paramIn = track->GetInnerParam();
524
525     // long track ..
526     if (!paramIn) continue;
527     if (!paramOut) continue;
528
529     // not a kink
530     if (track->GetKinkIndex(0) > 0) continue; 
531
532     Double_t extZ = GetExtZ(paramIn);
533     if (TMath::Abs(extZ) > 320) continue; // acceptance cut
534
535     // .. in the acceptance
536     Int_t sector = GetSector(paramOut->GetAlpha());
537     Int_t stack = GetStack(paramOut);
538
539     UInt_t u = 1;
540     UInt_t status = track->GetStatus();
541     for(Int_t bit=0; bit<32; bit++) 
542       if (u<<bit & status) FillESDsData(2,bit);
543
544     const Int_t knbits = 6; 
545     Int_t bit[6] = {0,0,0,0,0,0};    
546     bit[0] = status & AliESDtrack::kTPCin;
547     bit[1] = status & AliESDtrack::kTPCout;
548     bit[2] = (status & AliESDtrack::kTPCout) && !(status & AliESDtrack::kTRDout);
549     bit[3] = status & AliESDtrack::kTRDout;
550     bit[4] = status & AliESDtrack::kTRDrefit;
551     bit[5] = (status & AliESDtrack::kTRDout) && !(status & AliESDtrack::kTRDrefit);
552
553     // transverse momentum
554     //const Double_t *val = paramOut->GetParameter(); // parameters at the Outer plane
555     Double_t pt = paramOut->Pt(); //1./TMath::Abs(val[4]);
556
557     for(Int_t b=0; b<knbits; b++) {
558       if (bit[b]) {
559         FillESDsData(2*b+3,pt); 
560         FillESDsData(2*b+4,extZ);
561       }
562     }
563
564     // clusters
565     for(Int_t b=0; b<3; b++) 
566       if (bit[3+b]) FillESDsData(b+15,track->GetTRDncls0());
567
568     // refitted only
569     if (!bit[4]) continue;
570
571     //fQuality->Fill(track->GetTRDQuality());
572     //fBudget->Fill(track->GetTRDBudget());
573     //fSignal->Fill(track->GetTRDsignal());
574
575     FillESDsData(1,sector);
576     FillESDsData(18,track->GetP(), track->GetTRDsignal());
577
578     FillESDsData(27,track->GetTRDQuality());
579     FillESDsData(28,track->GetTRDBudget());
580     FillESDsData(29,track->GetTRDchi2());
581     FillESDsData(30,track->GetTRDTimBin(0));
582     FillESDsData(31,track->GetTRDntrackletsPID());
583     
584     
585     // dedx
586     for(Int_t k=0; k<4; ++k) {
587       Double_t dedx = 0;
588       for(Int_t j=0; j<6; j++) {
589         dedx += track->GetTRDslice(j, k-1);
590       }
591       FillESDsData(41+k,paramOut->GetP(), dedx/6.);
592     }
593
594     // probabilities
595     if (status & AliESDtrack::kTRDpid) {
596       for(Int_t k=0; k<AliPID::kSPECIES; ++k) 
597         FillESDsData(36+k,track->GetTRDpid(k));
598     }
599
600     // probabilities uniformity
601     if (track->GetTRDntrackletsPID() < 6) continue;
602     FillESDsData(35,paramOut->GetZ()/paramOut->GetX());
603     
604     Int_t idx = 5 * sector + stack;
605     FillESDsData(32,idx); // all tracks
606     if (track->GetTRDpid(AliPID::kElectron) > 0.9) 
607       FillESDsData(33,idx); // electrons only
608
609     
610
611     /*
612     hist[27] = new TH1F("qaTRD_esd_quality", ";quality", 120, 0, 12);
613     hist[28] = new TH1F("qaTRD_esd_budget", ";NN", 110, -1000, 100);
614     hist[29] = new TH1F("qaTRD_esd_chi2", ";chi2", 300, 0, 100);
615     hist[30] = new TH1F("qaTRD_esd_timeBin", 7, -0.5, 6.5);
616     hist[31] = new TH1F("qaTRD_esd_pidQuality", 7, -0.5, 6.5);
617     */
618
619     /*
620     // PID only
621     if (status & AliESDtrack::kTRDpid) {
622
623       for(Int_t l=0; l<6; l++) fTime->Fill(track->GetTRDTimBin(l));
624
625       // fill pid histograms
626       Double_t trdr0 = 0; //, tpcr0 = 0;
627       Int_t trdBestPid = 5; //, tpcBestPid = 5;  // charged
628       const Double_t kminPidValue = 0.9;
629
630       //Double_t pp[5];
631       //track->GetTPCpid(pp); // ESD inconsequence
632
633       for(Int_t pid=0; pid<5; pid++) {
634         
635         trdr0 += track->GetTRDpid(pid);
636         //tpcr0 += pp[pid];
637         
638         fTrdPID[pid]->Fill(track->GetTRDpid(pid));
639         //fTpcPID[pid]->Fill(pp[pid]);
640         
641         if (track->GetTRDpid(pid) > kminPidValue) trdBestPid = pid;
642         //if (pp[pid] > kminPidValue) tpcBestPid = pid;
643       }
644
645       fTrdPID[5]->Fill(trdr0); // check unitarity
646       fTrdSigMomPID[trdBestPid]->Fill(track->GetP(), track->GetTRDsignal());
647
648       //fTpcPID[5]->Fill(tpcr0); // check unitarity
649       //fTpcSigMomPID[tpcBestPid]->Fill(track->GetP(), track->GetTPCsignal());
650     }
651     */
652
653   }
654   //
655   IncEvCountCycleESDs();
656   IncEvCountTotalESDs();
657   //
658 }
659
660 //______________________________________________________________________________
661 Int_t AliTRDQADataMakerRec::GetSector(Double_t alpha) const 
662 {
663   //
664   // Gets the sector number 
665   //
666
667   Double_t size = TMath::DegToRad() * 20.; // shall use TRDgeo
668   if (alpha < 0) alpha += 2*TMath::Pi();
669   Int_t sector = (Int_t)(alpha/size);
670   return sector;
671
672 }
673 //______________________________________________________________________________
674
675 Int_t AliTRDQADataMakerRec::GetStack(const AliExternalTrackParam *paramOut) const
676 {
677   //
678   // calculates the stack the track is in
679   //
680   
681   const Double_t L = -0.9;
682   const Double_t W = (2*L)/5;
683
684   Double_t tan = paramOut->GetZ() / paramOut->GetX();
685   Double_t pos = (tan - L) / W;
686   return (Int_t) pos;
687 }
688
689 //______________________________________________________________________________
690 Double_t AliTRDQADataMakerRec::GetExtZ(const AliExternalTrackParam *in) const 
691 {
692   //
693   // Returns the Z position at the entry to TRD
694   // using parameters from the TPC in
695   //
696
697   const Double_t kX0 = 300;
698
699   Double_t x = in->GetX();
700   const Double_t *par = in->GetParameter();
701   Double_t theta = par[3];
702   Double_t z = in->GetZ();
703
704   Double_t zz = z + (kX0-x) * TMath::Tan(theta);
705   return zz;
706
707 }
708
709 //____________________________________________________________________________
710 void AliTRDQADataMakerRec::MakeRaws(AliRawReader* rawReader)
711 {
712   //
713   // Makes QA data from raw data
714   //
715
716   AliInfo("Making QA for Raws");
717
718   // take histograms (RS: arrays of histos)
719   static TObjArray mnADC,mnCls,mClsDet,mADC,mCls,mClsTb,mClsDetAmp;
720   GetMatchingRawsData(0,&mnADC);
721   GetMatchingRawsData(1,&mnCls); 
722   GetMatchingRawsData(2,&mClsDet);
723   GetMatchingRawsData(3,&mADC);
724   GetMatchingRawsData(4,&mCls);
725   GetMatchingRawsData(5,&mClsTb);
726   GetMatchingRawsData(6,&mClsDetAmp);
727
728   const Int_t baseline = 10;
729
730   // configure the reader
731   rawReader->Reset();
732   rawReader->SelectEquipment(0, 1024, 1041);
733   rawReader->Select("TRD");
734
735   AliTRDrawStream *data = new AliTRDrawStream(rawReader);
736
737   // build data manager  
738   AliTRDdigitsManager *digitsManager;
739   digitsManager = new AliTRDdigitsManager(kTRUE);
740   digitsManager->CreateArrays();
741
742   // error container 
743   const UShort_t kErrorChmb = 1411;
744   UShort_t **fErrorContainer = new UShort_t *[2];
745   fErrorContainer[0] = new UShort_t[kErrorChmb];
746   fErrorContainer[1] = new UShort_t[kErrorChmb];
747   
748   Int_t det = 0;
749   Int_t row, col;
750
751   while ((det = data->NextChamber(digitsManager, NULL, fErrorContainer)) >= 0){
752     
753     //printf("DET = %d\n", det);
754
755     AliTRDSignalIndex* indexes = digitsManager->GetIndexes(det);
756     if (indexes->HasEntry()) {
757       
758       AliTRDarrayADC *digits = (AliTRDarrayADC*) digitsManager->GetDigits(det);
759
760       while(indexes->NextRCIndex(row, col))  {
761
762         for (int ih=mnADC.GetEntriesFast();ih--;) ((TH2*)mnADC.UncheckedAt(ih))->Fill(det/30, det%30);
763         
764         for(Int_t tb = 0; tb < digits->GetNtime(); tb++) {
765           Int_t value = digits->GetData(row, col, tb);
766           for (int ih=mADC.GetEntriesFast();ih--;) ((TH1*)mADC.UncheckedAt(ih))->Fill(value);
767           
768           // simple clusterizer
769           if (col < 1 || col > digits->GetNcol()-2) continue;
770           if (tb < 1 || tb > digits->GetNtime()-2) continue;
771           
772           value -= baseline;
773
774           Int_t valueL = digits->GetData(row, col-1, tb) - baseline;
775           if (valueL >= value) continue;
776           
777           Int_t valueR = digits->GetData(row, col+1, tb) - baseline;
778           if (valueR >= value) continue;
779           
780           Int_t valueUp = digits->GetData(row, col-1, tb+1) + 
781             digits->GetData(row, col, tb+1) + digits->GetData(row, col+1, tb+1) - 3 * baseline;
782           if (valueUp < 10) continue;
783           
784           Int_t valueDown = digits->GetData(row, col-1, tb-1) + 
785             digits->GetData(row, col, tb-1) + digits->GetData(row, col+1, tb-1) - 3 * baseline;
786           if (valueDown < 10) continue;
787           
788           Int_t valueTot = value + valueL + valueR;
789           if (valueTot < 0) continue;
790
791           for (int ih=mCls.GetEntriesFast();ih--;) ((TH1*)mCls.UncheckedAt(ih))->Fill(valueTot);
792           for (int ih=mClsTb.GetEntriesFast();ih--;) ((TH2*)mClsTb.UncheckedAt(ih))->Fill(tb, valueTot);
793           for (int ih=mClsDetAmp.GetEntriesFast();ih--;) ((TH2*)mClsDetAmp.UncheckedAt(ih))->Fill(det, valueTot);
794
795           if (valueTot < 200) {
796             for (int ih=mnCls.GetEntriesFast();ih--;) ((TH2*)mnCls.UncheckedAt(ih))->Fill(det/30, det%30); 
797             for (int ih=mClsDet.GetEntriesFast();ih--;) ((TH2*)mClsDet.UncheckedAt(ih))->Fill(det/30, det%30, valueTot); 
798           }
799         
800         }
801       }
802
803       digitsManager->ClearArrays(det); // do we need this if object will be deleted ??
804     }    
805   }
806   
807   delete [] fErrorContainer[0];
808   delete [] fErrorContainer[1];
809   delete [] fErrorContainer;
810   fErrorContainer = NULL;
811   
812   delete digitsManager;  
813   delete data;
814   //
815   IncEvCountCycleRaws();
816   IncEvCountTotalRaws();
817   //
818 }
819
820 //____________________________________________________________________________
821 void AliTRDQADataMakerRec::MakeRecPoints(TTree * clustersTree)
822 {
823   //  
824   // Makes data from RecPoints
825   // 
826   
827   //  Info("MakeRecPoints", "making");
828  
829   Int_t nsize = Int_t(clustersTree->GetTotBytes() / (sizeof(AliTRDcluster))); 
830   TObjArray *clusterArray = new TObjArray(nsize+1000); 
831
832   TBranch *branch = clustersTree->GetBranch("TRDcluster");
833   if (!branch) {
834     AliError("Can't get the branch !");
835     return;
836   }
837   branch->SetAddress(&clusterArray); 
838
839   // Loop through all entries in the tree
840   Int_t nEntries   = (Int_t)TMath::Ceil( clustersTree->GetEntries() );
841   Int_t nbytes     = 0;
842   AliTRDcluster *c = 0;
843   Int_t nDet[540];
844   for (Int_t i=0; i<540; i++) nDet[i] = 0;
845
846   Int_t nCls = 0;
847   
848   //printf("nEntries = %d\n", nEntries);
849   //nEntries++;
850   
851   /*
852   // select the event 
853   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
854
855     // Import the tree
856     nbytes += clustersTree->GetEvent(iEntry);  
857     Int_t nCluster = clusterArray->GetEntries();  
858     
859     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
860       c = (AliTRDcluster *) clusterArray->At(iCluster);
861       nCls++;
862     }
863   }
864
865   if (nCls < 100) {
866     delete clusterArray;
867     return;
868   }
869   */
870
871   /////
872   TObjArray *hists3D = GetMatchingRecPointsData(10);       //RS no alias for 3d histo filling, to directly
873   //
874   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
875
876     // Import the tree
877     nbytes += clustersTree->GetEvent(iEntry);  
878
879     // Get the number of points in the detector
880     Int_t nCluster = clusterArray->GetEntries();  
881
882   
883     // Loop through all TRD digits
884     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
885
886       nCls++;
887       c = (AliTRDcluster *) clusterArray->At(iCluster);
888
889       Int_t iDet = c->GetDetector();
890       Int_t nPads = c->GetNPads();
891
892       nDet[iDet]++;
893       FillRecPointsData(0,iDet);
894       FillRecPointsData(86,c->GetQ());
895       FillRecPointsData(1,iDet, c->GetQ());
896       FillRecPointsData(2,nPads);
897       if (nPads < 6)
898         FillRecPointsData(1+c->GetNPads(),c->GetCenter());
899       
900       if (nPads < 10)
901         FillRecPointsData(88+nPads-2,c->GetQ());
902       else FillRecPointsData(96,c->GetQ());
903
904       //if (c->GetPadTime() < 5)
905       FillRecPointsData(7,c->GetPadRow(), c->GetPadCol());
906       FillRecPointsData(8,c->GetPadTime());
907       //
908       if (hists3D) {
909         for (int ih=hists3D->GetEntriesFast();ih--;) {
910           TH3F * ahist = dynamic_cast<TH3F*>(hists3D->At(ih));
911           if (ahist) ahist->Fill(iDet, c->GetPadTime(), c->GetQ());
912         }
913       }
914       
915       Int_t iSM = iDet / 30;
916       FillRecPointsData(50+iSM,c->GetPadTime());
917       FillRecPointsData(68+iSM,c->GetPadTime(), c->GetQ());
918       
919       // total charge sm / det / timeBin
920       // FillRecPointsData(14+iSM,iDet-iSM*30, c->GetPadTime(), c->GetQ());
921
922
923       // PRF for 2pad
924       //if (c->GetNPads() == 2) {
925       Short_t *sig = c->GetSignals();
926       Double_t frac = -10;
927
928       if (sig[0] == 0 && sig[1] == 0 && sig[2] == 0 && sig[5] == 0 && sig[6] == 0) 
929         frac = 1. * sig[4] / (sig[3] + sig[4]);
930
931       if (sig[0] == 0 && sig[1] == 0 && sig[4] == 0 && sig[5] == 0 && sig[6] == 0)
932         frac = -1. * sig[2] / (sig[2] + sig[3]);
933
934       if (frac > -10) FillRecPointsData(11,c->GetCenter(), frac);
935         
936       //}
937     }
938     clusterArray->Delete();
939   }
940
941   /*
942   for(Int_t i=0; i<540; i++) 
943     if (nDet[i] > 0) FillRecPointsData(9,nDet[i]);
944   */
945   FillRecPointsData(9,nCls);
946   
947   delete clusterArray;
948   //
949  IncEvCountCycleRecPoints();
950  IncEvCountTotalRecPoints();
951   //
952 }
953
954 //____________________________________________________________________________ 
955 void AliTRDQADataMakerRec::StartOfDetectorCycle()
956 {
957   //
958   // Detector specific actions at start of cycle
959   //
960
961 }
962
963 //__________________________________________________________________________
964 Int_t AliTRDQADataMakerRec::CheckPointer(TObject *obj, const char *name) 
965 {
966   //
967   // Checks initialization of pointers
968   //
969
970   if (!obj) AliWarning(Form("null pointer: %s", name));
971   return !!obj;
972 }
973 //__________________________________________________________________________
974 void AliTRDQADataMakerRec::BuildRatio(TH1 *ratio, TH1 *histN, TH1* histD) {
975   //
976   // Calculate the ratio of two histograms 
977   // error are calculated assuming the histos have the same counts
978   //
979
980   // calclate
981
982   Int_t nbins = histN->GetXaxis()->GetNbins();
983   for(Int_t i=1; i<nbins+2; i++) {
984     
985     Double_t valueN = histN->GetBinContent(i);
986     Double_t valueD = histD->GetBinContent(i);
987     
988     if (valueD < 1) {
989       ratio->SetBinContent(i, 0);
990       ratio->SetBinError(i, 0);
991       continue;
992     }
993
994     Double_t eps = (valueN < valueD-valueN)? valueN : valueD-valueN;
995     
996     ratio->SetBinContent(i, valueN/valueD);
997     ratio->SetBinError(i, TMath::Sqrt(eps)/valueD);
998   }
999
1000   // style
1001   ratio->SetMinimum(-0.1);
1002   ratio->SetMaximum(1.1);
1003   ratio->SetMarkerStyle(20);
1004 }
1005 //__________________________________________________________________________
1006
1007 Int_t AliTRDQADataMakerRec::FillBits(TH1F *hist, Int_t code, Int_t offset) {
1008
1009   Int_t nb = 0;
1010   UInt_t test = 1;
1011   for(Int_t i=0; i<8; i++) {
1012     if (code & test) {
1013       hist->Fill(i+offset);
1014       nb++;
1015     }
1016     test *= 2;       
1017   }
1018   
1019   return nb;
1020 }
1021
1022 //__________________________________________________________________________