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