TPC/AliTPCQADataMakerRec.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, 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
  17 /* $Id: $ */
  18
  19 /*
  20   Based on AliPHOSQADataMaker
  21   Produces the data needed to calculate the quality assurance.
  22   All data must be mergeable objects.
  23   P. Christiansen, Lund, January 2008
  24 */
  25
  26 /*
  27   Implementation:
  28
  29   We have chosen to have the histograms as non-persistent meber to
  30   allow better debugging. In the copy constructor we then have to
  31   assign the pointers to the existing histograms in the copied
  32   list. This have been implemented but not tested.
  33
  34   For the QA of the RAW data we use the class, AliTPCdataQA, from the
  35   existing TPC Calibration framework (which is more advanced than the
  36   standard QA framework) and extract the histograms at the end. This
  37   has been tested with zero-suppressed data. The Analyse method of the
  38   AliTPCdataQA class is called in the method, EndOfDetectorCycle, and
  39   there also: 1d histogram(s) are projected and added to the QA list.
  40 */
  41
  42 /*
  43   TODO:
  44   Sumw2 for RAW histogram(s)?
  45   RecPoints and ESD could have many more histograms
  46 */
  47
  48 #include "AliTPCQADataMakerRec.h"
  49
  50 // --- ROOT system ---
  51 #include <TClonesArray.h>
  52 #include <TString.h>
  53 #include <TSystem.h>
  54
  55 // --- Standard library ---
  56
  57 // --- AliRoot header files ---
  58 #include "AliQAChecker.h"
  59 #include "AliESDEvent.h"
  60 #include "AliESDtrack.h"
  61 #include "AliLog.h"
  62 #include "AliTPCCalPad.h"
  63 #include "AliTPCCalROC.h"
  64 #include "AliTPCClustersRow.h"
  65 #include "AliTPCclusterMI.h"
  66 #include "AliSimDigits.h"
  67
  68 ClassImp(AliTPCQADataMakerRec)
  69
  70 //____________________________________________________________________________
  71 AliTPCQADataMakerRec::AliTPCQADataMakerRec() :
  72   AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kTPC),
  73                     "TPC Rec Quality Assurance Data Maker"),
  74   fTPCdataQA(NULL)
  75 {
  76   // ctor
  77   fTPCdataQA = new AliTPCdataQA*[AliRecoParam::kNSpecies] ;
  78   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++)
  79     fTPCdataQA[specie] = NULL ;
  80
  81   for(Int_t i = 0; i < 6; i++)
  82     fMapping[i] = 0;
  83 }
  84
  85 //____________________________________________________________________________
  86 AliTPCQADataMakerRec::AliTPCQADataMakerRec(const AliTPCQADataMakerRec& qadm) :
  87   AliQADataMakerRec(),
  88   fTPCdataQA(NULL)
  89 {
  90   //copy ctor
  91   // Does not copy the calibration object, instead InitRaws have to be
  92   // called again
  93   SetName((const char*)qadm.GetName()) ;
  94   SetTitle((const char*)qadm.GetTitle());
  95
  96   fTPCdataQA = new AliTPCdataQA*[AliRecoParam::kNSpecies] ;
  97   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++)
  98     fTPCdataQA[specie] = NULL ;
  99
 100   for(Int_t i = 0; i < 6; i++)
 101     fMapping[i] = 0;
 102
 103   //
 104   // Associate class histogram objects to the copies in the list
 105   // Could also be done with the indexes
 106   //
 107
 108 }
 109
 110 //__________________________________________________________________
 111 AliTPCQADataMakerRec& AliTPCQADataMakerRec::operator = (const AliTPCQADataMakerRec& qadm )
 112 {
 113   // Equal operator.
 114   this->~AliTPCQADataMakerRec();
 115   new(this) AliTPCQADataMakerRec(qadm);
 116   return *this;
 117 }
 118
 119 //__________________________________________________________________
 120 AliTPCQADataMakerRec::~AliTPCQADataMakerRec()
 121 {
 122   // Destructor
 123   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++)
 124     if ( fTPCdataQA[specie] != NULL )
 125       delete fTPCdataQA[specie] ;
 126   delete[] fTPCdataQA;
 127
 128   for(Int_t i = 0; i < 6; i++)
 129     delete fMapping[i];
 130 }
 131
 132 //____________________________________________________________________________
 133 void AliTPCQADataMakerRec::EndOfDetectorCycle(AliQAv1::TASKINDEX_t task, TObjArray ** list)
 134 {
 135   //Detector specific actions at end of cycle
 136
 137   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
 138     if ( !AliQAv1::Instance()->IsEventSpecieSet(specie) )
 139       continue ;
 140     if(fTPCdataQA[specie] != NULL) { // do the final step of the QA for Raw data
 141
 142       fTPCdataQA[specie]->Analyse(); // 31/1-08 Analyse is now protected against
 143                            //         RAW data files with no TPC data
 144
 145       SetEventSpecie(AliRecoParam::ConvertIndex(specie)) ;
 146       TH1F * histRawsOccupancy                 = (TH1F*)GetRawsData(kOccupancy) ;
 147       TH1F * histRawsOccupancyVsSector         = (TH1F*)GetRawsData(kOccupancyVsSector) ;
 148       TH1F * histRawsNClustersPerEventVsSector = (TH1F*)GetRawsData(kNClustersPerEventVsSector) ;
 149       TH1F * histRawsQVsSector                 = (TH1F*)GetRawsData(kQVsSector) ;
 150       TH1F * histRawsQmaxVsSector              = (TH1F*)GetRawsData(kQmaxVsSector) ;
 151       if ( !histRawsOccupancy ||
 152           !histRawsOccupancyVsSector ||
 153           !histRawsNClustersPerEventVsSector ||
 154           !histRawsQVsSector ||
 155           !histRawsQmaxVsSector) {
 156         AliError("Something very wrong here, corrupted memory ?????. Please check\n") ;
 157         continue ;
 158       }
 159
 160       //Add2RawsList(fTPCdataQA, 0);
 161       // get the histograms and add them to the output
 162       // 31/8-08 Histogram is only added if the Calibration class
 163       //         receives TPC data
 164       const Int_t eventCounter = fTPCdataQA[specie]->GetEventCounter();
 165       if(eventCounter>0) { // some TPC data has been processed
 166
 167         // Reset histograms and refill them
 168         histRawsOccupancy->Reset();
 169         histRawsOccupancyVsSector->Reset();
 170         histRawsNClustersPerEventVsSector->Reset();
 171         histRawsQVsSector->Reset();
 172         histRawsQmaxVsSector->Reset();
 173
 174         TH1F* hNorm72 = new TH1F("hNorm72", "histogram to normalize 72 sectors",
 175                                  72, 0, 72);
 176         hNorm72->Sumw2();
 177         TH1F* hNorm108 = new TH1F("hNorm108", "histogram to normalize 108 sectors (medium and long pads are split up)",
 178                                   108, 0, 108);
 179         hNorm108->Sumw2();
 180
 181         for (Int_t iSec = 0; iSec < 72; iSec++) {
 182
 183           AliTPCCalROC* occupancyROC =
 184           fTPCdataQA[specie]->GetNoThreshold()->GetCalROC(iSec);
 185           AliTPCCalROC* nclusterROC =
 186           fTPCdataQA[specie]->GetNLocalMaxima()->GetCalROC(iSec);
 187           AliTPCCalROC* qROC =
 188           fTPCdataQA[specie]->GetMeanCharge()->GetCalROC(iSec);
 189           AliTPCCalROC* qmaxROC =
 190           fTPCdataQA[specie]->GetMaxCharge()->GetCalROC(iSec);
 191
 192           const Int_t nRows = occupancyROC->GetNrows();
 193           for (Int_t iRow = 0; iRow < nRows; iRow++) {
 194
 195             Int_t helpSector = iSec;
 196             if(iRow>=64)
 197               helpSector += 36; // OROC (long pads)
 198
 199             const Int_t nPads = occupancyROC->GetNPads(iRow);
 200             for (Int_t iPad = 0; iPad < nPads; iPad++) {
 201
 202               histRawsOccupancy->Fill(occupancyROC->GetValue(iRow, iPad));
 203               hNorm72->Fill(iSec);
 204               histRawsOccupancyVsSector
 205               ->Fill(iSec, occupancyROC->GetValue(iRow, iPad));
 206
 207               const Int_t nClusters = TMath::Nint(nclusterROC->GetValue(iRow, iPad));
 208
 209               if(nClusters>0) {
 210
 211                 histRawsNClustersPerEventVsSector->Fill(iSec, nClusters);
 212                 hNorm108->Fill(helpSector, nClusters);
 213                 histRawsQVsSector->Fill(helpSector,
 214                                        nClusters*qROC->GetValue(iRow, iPad));
 215                 histRawsQmaxVsSector->Fill(helpSector,
 216                                           nClusters*qmaxROC->GetValue(iRow, iPad));
 217               }
 218             }
 219           }
 220         } // end loop over sectors
 221
 222         // Normalize histograms
 223         histRawsOccupancyVsSector->Divide(hNorm72);
 224         histRawsNClustersPerEventVsSector->Scale(1.0/Float_t(eventCounter));
 225         histRawsQVsSector->Divide(hNorm108);
 226         histRawsQmaxVsSector->Divide(hNorm108);
 227         delete hNorm72;
 228         delete hNorm108;
 229       }
 230     }
 231   }
 232   AliQAChecker::Instance()->Run(AliQAv1::kTPC, task, list) ;
 233 }
 234
 235 //____________________________________________________________________________
 236 void AliTPCQADataMakerRec::InitESDs()
 237 {
 238   //create ESDs histograms in ESDs subdir
 239   const Bool_t expert   = kTRUE ;
 240   const Bool_t image    = kTRUE ;
 241
 242   TH1F * histESDclusters =
 243     new TH1F("hESDclusters", "N TPC clusters per track; N clusters; Counts",
 244              160, 0, 160);
 245   histESDclusters->Sumw2();
 246   Add2ESDsList(histESDclusters, KClusters, !expert, image);
 247
 248   TH1F * histESDratio =
 249     new TH1F("hESDratio", "Ratio: TPC clusters / findable; Ratio: cluster/findable; Counts",
 250              100, 0, 1);
 251   histESDratio->Sumw2();
 252   Add2ESDsList(histESDratio, kRatio, !expert, image);
 253
 254   TH1F * histESDpt =
 255     new TH1F("hESDpt", "P_{T} distribution; p_{T} [GeV/c]; Counts",
 256              50, 0, 5);
 257   histESDpt->Sumw2();
 258   Add2ESDsList(histESDpt, kPt, !expert, image);
 259 }
 260
 261 //____________________________________________________________________________
 262 void AliTPCQADataMakerRec::InitRaws()
 263 {
 264   //
 265   // Adding the raw
 266   //
 267
 268   // Modified: 7/7 - 2008
 269   // Laurent Aphecetche pointed out that the mapping was read from file
 270   // for each event, so now we read in the map here and set if for
 271   // the raw data qa
 272   const Bool_t expert   = kTRUE ;
 273   const Bool_t saveCorr = kTRUE ;
 274   const Bool_t image    = kTRUE ;
 275
 276   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
 277     fTPCdataQA[specie] = new AliTPCdataQA(AliRecoParam::ConvertIndex(specie));
 278     LoadMaps(); // Load Altro maps
 279     fTPCdataQA[specie]->SetAltroMapping(fMapping); // set Altro mapping
 280     fTPCdataQA[specie]->SetRangeTime(100, 920); // set time bin interval
 281 //    Add2RawsList(fTPCdataQA, kTPCdataQ, !expert, image, !saveCorrA); // This is used by the AMORE monitoring <------- THIS WILL FAIL (YS)
 282   }
 283
 284   TH1F * histRawsOccupancy =
 285     new TH1F("hRawsOccupancy", "Occupancy (all pads); Occupancy; Counts",
 286              100, 0, 1);
 287   histRawsOccupancy->Sumw2();
 288   Add2RawsList(histRawsOccupancy, kOccupancy, !expert, image, !saveCorr);
 289
 290   TH1F * histRawsOccupancyVsSector =
 291     new TH1F("hRawsOccupancyVsSector", "Occupancy vs sector; Sector; Occupancy",
 292              72, 0, 72);
 293   histRawsOccupancyVsSector->Sumw2();
 294   Add2RawsList(histRawsOccupancyVsSector, kOccupancyVsSector, !expert, image, !saveCorr);
 295
 296   TH1F * histRawsNClustersPerEventVsSector =
 297     new TH1F("hRawsNClustersPerEventVsSector", "Nclusters per event vs sector; Sector; Nclusters per event",
 298              72, 0, 72);
 299   histRawsNClustersPerEventVsSector->Sumw2();
 300   Add2RawsList(histRawsNClustersPerEventVsSector, kNClustersPerEventVsSector, !expert, image, !saveCorr);
 301
 302   TH1F * histRawsQVsSector =
 303     new TH1F("hRawsQVsSector", "<Q> vs sector (OROC med: 36-71, long: 72-107); Sector; <Q>",
 304              108, 0, 108);
 305   histRawsQVsSector->Sumw2();
 306   Add2RawsList(histRawsQVsSector, kQVsSector, !expert, image, !saveCorr);
 307
 308   TH1F * histRawsQmaxVsSector =
 309     new TH1F("hRawsQmaxVsSector", "<Qmax> vs sector (OROC med: 36-71, long: 72-107); Sector; <Qmax>",
 310              108, 0, 108);
 311   histRawsQmaxVsSector->Sumw2();
 312   Add2RawsList(histRawsQmaxVsSector, kQmaxVsSector, !expert, image, !saveCorr);
 313 }
 314
 315 //____________________________________________________________________________
 316 void AliTPCQADataMakerRec::InitDigits()
 317 {
 318   const Bool_t expert   = kTRUE ;
 319   const Bool_t image    = kTRUE ;
 320   TH1F * histDigitsADC =
 321     new TH1F("hDigitsADC", "Digit ADC distribution; ADC; Counts",
 322              1000, 0, 1000);
 323   histDigitsADC->Sumw2();
 324   Add2DigitsList(histDigitsADC, kDigitsADC, !expert, image);
 325 }
 326
 327 //____________________________________________________________________________
 328 void AliTPCQADataMakerRec::InitRecPoints()
 329 {
 330   const Bool_t expert   = kTRUE ;
 331   const Bool_t image    = kTRUE ;
 332
 333   TH1F * histRecPointsQmaxShort =
 334     new TH1F("hRecPointsQmaxShort", "Qmax distrbution (short pads); Qmax; Counts",
 335              100, 0, 300);
 336   histRecPointsQmaxShort->Sumw2();
 337   Add2RecPointsList(histRecPointsQmaxShort, kQmaxShort, !expert, image);
 338
 339   TH1F * histRecPointsQmaxMedium =
 340     new TH1F("hRecPointsQmaxMedium", "Qmax distrbution (medium pads); Qmax; Counts",
 341              100, 0, 300);
 342   histRecPointsQmaxMedium->Sumw2();
 343   Add2RecPointsList(histRecPointsQmaxMedium, kQmaxMedium, !expert, image);
 344
 345   TH1F * histRecPointsQmaxLong =
 346     new TH1F("hRecPointsQmaxLong", "Qmax distrbution (long pads); Qmax; Counts",
 347              100, 0, 300);
 348   histRecPointsQmaxLong->Sumw2();
 349   Add2RecPointsList(histRecPointsQmaxLong, kQmaxLong, !expert, image);
 350
 351   TH1F * histRecPointsQShort =
 352     new TH1F("hRecPointsQShort", "Q distrbution (short pads); Q; Counts",
 353              100, 0, 2000);
 354   histRecPointsQShort->Sumw2();
 355   Add2RecPointsList(histRecPointsQShort, kQShort, !expert, image);
 356
 357   TH1F * histRecPointsQMedium =
 358     new TH1F("hRecPointsQMedium", "Q distrbution (medium pads); Q; Counts",
 359              100, 0, 2000);
 360   histRecPointsQMedium->Sumw2();
 361   Add2RecPointsList(histRecPointsQMedium, kQMedium, !expert, image);
 362
 363   TH1F * histRecPointsQLong =
 364     new TH1F("hRecPointsQLong", "Q distrbution (long pads); Q; Counts",
 365              100, 0, 2000);
 366   histRecPointsQLong->Sumw2();
 367   Add2RecPointsList(histRecPointsQLong, kQLong, !expert, image);
 368
 369   TH1F * histRecPointsRow =
 370     new TH1F("hRecPointsRow", "Clusters per row; Row; Counts",
 371              159, 0, 159);
 372   histRecPointsRow->Sumw2();
 373   Add2RecPointsList(histRecPointsRow, kRow, !expert, image);
 374 }
 375
 376 //____________________________________________________________________________
 377 void AliTPCQADataMakerRec::MakeESDs(AliESDEvent * esd)
 378 {
 379   // make QA data from ESDs
 380
 381   const Int_t nESDTracks = esd->GetNumberOfTracks();
 382   Int_t nTPCtracks = 0;
 383   for(Int_t i = 0; i < nESDTracks; i++) {
 384
 385     AliESDtrack * track = esd->GetTrack(i);
 386
 387     if ((track->GetStatus() & AliESDtrack::kTPCrefit)==0)
 388       continue;
 389
 390     nTPCtracks++;
 391
 392     Int_t nTPCclusters         = track->GetTPCNcls();
 393     Int_t nTPCclustersFindable = track->GetTPCNclsF();
 394     if ( nTPCclustersFindable<=0) continue;
 395     GetESDsData(KClusters)->Fill(nTPCclusters);
 396     GetESDsData(kRatio)->Fill(Float_t(nTPCclusters)/Float_t(nTPCclustersFindable));
 397     GetESDsData(kPt)->Fill(track->Pt());
 398   }
 399 }
 400
 401 //____________________________________________________________________________
 402 void AliTPCQADataMakerRec::MakeRaws(AliRawReader* rawReader)
 403 {
 404   //
 405   // To make QA for the RAW data we use the TPC Calibration framework
 406   // to handle the data and then in the end extract the data
 407   //
 408
 409   GetRawsData(0); // dummy call to init raw data
 410   rawReader->Reset() ;
 411   if (! fTPCdataQA[AliRecoParam::AConvert(fEventSpecie)] ) {
 412     AliError("Something unexpected here!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!") ;
 413   } else {
 414     fTPCdataQA[AliRecoParam::AConvert(fEventSpecie)]->ProcessEvent(rawReader);
 415   }
 416 }
 417
 418 //____________________________________________________________________________
 419 void AliTPCQADataMakerRec::MakeDigits(TTree* digitTree)
 420 {
 421
 422   TBranch* branch = digitTree->GetBranch("Segment");
 423   AliSimDigits* digArray = 0;
 424   branch->SetAddress(&digArray);
 425
 426   Int_t nEntries = Int_t(digitTree->GetEntries());
 427
 428   for (Int_t n = 0; n < nEntries; n++) {
 429
 430     digitTree->GetEvent(n);
 431
 432     if (digArray->First())
 433       do {
 434         Float_t dig = digArray->CurrentDigit();
 435
 436         GetDigitsData(kDigitsADC)->Fill(dig);
 437       } while (digArray->Next());
 438   }
 439 }
 440
 441 //____________________________________________________________________________
 442 void AliTPCQADataMakerRec::MakeRecPoints(TTree* recTree)
 443 {
 444
 445   AliTPCClustersRow *clrow = new AliTPCClustersRow();
 446   clrow->SetClass("AliTPCclusterMI");
 447   clrow->SetArray(0);
 448   clrow->GetArray()->ExpandCreateFast(10000);
 449
 450   TBranch* branch = recTree->GetBranch("Segment");
 451   branch->SetAddress(&clrow);
 452
 453   const Int_t nEntries = Int_t(recTree->GetEntries());
 454   for (Int_t i = 0; i < nEntries; i++) {
 455
 456     branch->GetEntry(i);
 457
 458     const Int_t nClusters = clrow->GetArray()->GetEntriesFast();
 459     for (Int_t icl=0; icl < nClusters; icl++){
 460
 461       AliTPCclusterMI* cluster =
 462         (AliTPCclusterMI*)clrow->GetArray()->At(icl);
 463
 464       Float_t Qmax = cluster->GetMax();
 465       Float_t Q    = cluster->GetQ();
 466       Int_t   row  = cluster->GetRow();
 467
 468       if(cluster->GetDetector()<36) { // IROC (short pads)
 469
 470         GetRecPointsData(kQmaxShort)->Fill(Qmax);
 471         GetRecPointsData(kQShort)->Fill(Q);
 472       } else { // OROC (medium and long pads)
 473         row += 63;
 474         if(cluster->GetRow()<64) { // medium pads
 475
 476           GetRecPointsData(kQmaxMedium)->Fill(Qmax);
 477           GetRecPointsData(kQMedium)->Fill(Q);
 478         } else { // long pads
 479
 480           GetRecPointsData(kQmaxLong)->Fill(Qmax);
 481           GetRecPointsData(kQLong)->Fill(Q);
 482         }
 483       }
 484
 485       GetRecPointsData(kRow)->Fill(row);
 486     } // end loop over clusters
 487   } // end loop over tree
 488
 489   delete clrow;
 490 }
 491
 492 //____________________________________________________________________________
 493 void AliTPCQADataMakerRec::LoadMaps()
 494 {
 495   TString path = gSystem->Getenv("ALICE_ROOT");
 496   path += "/TPC/mapping/Patch";
 497
 498   for(Int_t i = 0; i < 6; i++) {
 499     TString path2 = path;
 500     path2 += i;
 501     path2 += ".data";
 502     fMapping[i] = new AliTPCAltroMapping(path2.Data());
 503   }
 504 }
 505