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 <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 //__________________________________________________________________________