TRD/AliTRDQADataMakerRec.cxx

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