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