TRD/qaRec/AliTRDcheckDetector.cxx

   1 #include <TFile.h>
   2 #include <TH1F.h>
   3 #include <TGraph.h>
   4 #include <TMath.h>
   5 #include <TMap.h>
   6 #include <TObjArray.h>
   7 #include <TObject.h>
   8 #include <TObjString.h>
   9 #include <TProfile.h>
  10 #include <TProfile2D.h>
  11 #include <TROOT.h>
  12
  13 #include "AliTRDcluster.h"
  14 #include "AliESDHeader.h"
  15 #include "AliESDRun.h"
  16 #include "AliTRDseedV1.h"
  17 #include "AliTRDtrackV1.h"
  18 #include "TTreeStream.h"
  19
  20 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  21 #include "AliTRDtrackInfo/AliTRDeventInfo.h"
  22 #include "AliTRDcheckDetector.h"
  23
  24 #include <cstdio>
  25 #include <iostream>
  26
  27 ////////////////////////////////////////////////////////////////////////////
  28 //                                                                        //
  29 //  Reconstruction QA                                                     //
  30 //                                                                        //
  31 //  Task doing basic checks for tracking and detector performance         //
  32 //                                                                        //
  33 //  Authors:                                                              //
  34 //    Anton Andronic <A.Andronic@gsi.de>                                  //
  35 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  36 //                                                                        //
  37 ////////////////////////////////////////////////////////////////////////////
  38
  39 //_______________________________________________________
  40 AliTRDcheckDetector::AliTRDcheckDetector():
  41   AliTRDrecoTask("DetChecker", "Basic Detector Checker")
  42   ,fEventInfo(0x0)
  43   ,fTriggerNames(0x0)
  44 {
  45   //
  46   // Default constructor
  47   //
  48   DefineInput(1,AliTRDeventInfo::Class());
  49 }
  50
  51 //_______________________________________________________
  52 AliTRDcheckDetector::~AliTRDcheckDetector(){
  53   //
  54   // Destructor
  55   //
  56   if(fEventInfo) delete fEventInfo;
  57   if(fTriggerNames) delete fTriggerNames;
  58 }
  59
  60 //_______________________________________________________
  61 void AliTRDcheckDetector::ConnectInputData(Option_t *opt){
  62         //
  63         // Connect the Input data with the task
  64         //
  65         AliTRDrecoTask::ConnectInputData(opt);
  66         fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
  67 }
  68
  69 //_______________________________________________________
  70 void AliTRDcheckDetector::CreateOutputObjects(){
  71   //
  72   // Create Output Objects
  73   //
  74   fContainer = new TObjArray(25);
  75   // Register Histograms
  76   fContainer->Add(new TH1F("hNtrks", "Number of Tracks per event", 100, 0, 100));
  77   fContainer->Add(new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100));
  78   fContainer->Add(new TH1F("hNcls", "Nr. of clusters per track", 181, -0.5, 180.5));
  79   fContainer->Add(new TH1F("hNtls", "Nr. tracklets per track", 7, -0.5, 6.5));
  80   fContainer->Add(new TH1F("hNclTls","Mean Number of clusters per tracklet", 31, -0.5, 30.5));
  81   fContainer->Add(new TH1F("hChi2", "Chi2", 200, 0, 20));
  82   fContainer->Add(new TH1F("hChi2n", "Norm. Chi2 (tracklets)", 50, 0, 5));
  83   fContainer->Add(new TH1F("hSM", "Track Counts in Supermodule", 18, -0.5, 17.5));
  84         // Detector signal on Detector-by-Detector basis
  85   fContainer->Add(new TProfile("hPHdetector", "Average PH", 31, -0.5, 30.5));
  86   fContainer->Add(new TH1F("hQclDetector", "Cluster charge", 200, 0, 1200));
  87   fContainer->Add(new TH1F("hQTdetector", "Total Charge Deposit", 6000, 0, 6000));
  88   fContainer->Add(new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100));
  89
  90   fTriggerNames = new TMap();
  91 }
  92
  93 //_______________________________________________________
  94 void AliTRDcheckDetector::Exec(Option_t *){
  95   //
  96   // Execution function
  97   // Filling TRD quality histos
  98   //
  99   if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
 100   Int_t nTracks = 0;            // Count the number of tracks per event
 101   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
 102   TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
 103   if(fDebugLevel > 6)printf("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data());
 104   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger))->Fill(triggermask);
 105   AliTRDtrackInfo *fTrackInfo = 0x0;
 106   AliTRDtrackV1 *fTRDtrack = 0x0;
 107   AliTRDseedV1 *fTracklet = 0x0;
 108   AliTRDcluster *fTRDcluster = 0x0;
 109   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
 110     fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
 111     if(!fTrackInfo || !(fTRDtrack = fTrackInfo->GetTRDtrack())) continue;
 112     Int_t nclusters = fTRDtrack->GetNumberOfClusters();
 113     Int_t ntracklets = fTRDtrack->GetNumberOfTracklets();
 114     Float_t chi2 = fTRDtrack->GetChi2();
 115     // Fill Histograms
 116     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist))->Fill(nclusters);
 117     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist))->Fill(ntracklets);
 118     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2))->Fill(chi2);
 119     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2Normalized))->Fill(chi2/static_cast<Float_t>(ntracklets));
 120     // now loop over single tracklets
 121     if(ntracklets > 2){
 122       Int_t sector = -1;
 123       for(Int_t ilayer = 0; ilayer < kNLayers; ilayer++){
 124         if(!(fTracklet = fTRDtrack->GetTracklet(ilayer)) || !fTracklet->IsOK()) continue;
 125         Float_t nClustersTracklet = fTracklet->GetN();
 126         dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist))->Fill(nClustersTracklet);
 127         if(nClustersTracklet){
 128           Float_t Qtot = 0;
 129           Int_t detector = -1;
 130           for(Int_t itb = 0; itb < kNTimebins; itb++){
 131             if(!(fTRDcluster = fTracklet->GetClusters(itb))) continue;
 132             Int_t localtime        = fTRDcluster->GetLocalTimeBin();
 133             Double_t clusterCharge = fTRDcluster->GetQ();
 134             detector               = fTRDcluster->GetDetector();
 135             sector                 = static_cast<Int_t>(detector/kNDetectorsSector);
 136             Double_t absolute_charge = TMath::Abs(clusterCharge);
 137             Qtot += absolute_charge;
 138             dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight))->Fill(localtime, absolute_charge);
 139             dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge))->Fill(absolute_charge);
 140             if(fDebugLevel > 2){
 141                 (*fDebugStream) << "PulseHeight"
 142                         << "Detector="  << detector
 143                         << "Sector="            << sector
 144                         << "Timebin="           << localtime
 145                         << "Charge="            << absolute_charge
 146                         << "\n";
 147             }
 148           }
 149           dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit))->Fill(Qtot);
 150           if(fDebugLevel > 3){
 151                 (*fDebugStream) << "ChargeDeposit"
 152                         << "Detector="  << detector
 153                         << "Sector="            << sector
 154                         << "QT="                                << Qtot
 155                         << "\n";
 156           }
 157         }
 158       }
 159       dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist))->Fill(sector);
 160     }
 161     nTracks++;
 162   }
 163   if(nTracks){
 164         dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks))->Fill(triggermask);
 165         dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist))->Fill(nTracks);
 166   }
 167   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
 168     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
 169     // also set the label for both histograms
 170     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
 171     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 172     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
 173     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 174   }
 175   PostData(0, fContainer);
 176 }
 177
 178 //_______________________________________________________
 179 void AliTRDcheckDetector::Terminate(Option_t *){
 180   //
 181   // Terminate function
 182   //
 183 }
 184
 185 //_______________________________________________________
 186 Bool_t AliTRDcheckDetector::PostProcess(){
 187         //
 188         // Do Postprocessing (for the moment set the number of Reference histograms)
 189         //
 190
 191         TH1 * histo = 0x0;
 192         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist));
 193         histo->GetXaxis()->SetTitle("Number of Tracks");
 194         histo->GetYaxis()->SetTitle("Events");
 195         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist));
 196         histo->GetXaxis()->SetTitle("Number of Clusters");
 197         histo->GetYaxis()->SetTitle("Entries");
 198         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist));
 199         histo->GetXaxis()->SetTitle("Number of Tracklets");
 200         histo->GetYaxis()->SetTitle("Entries");
 201         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist));
 202         histo->GetXaxis()->SetTitle("Number of Clusters");
 203         histo->GetYaxis()->SetTitle("Entries");
 204         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2));
 205         histo->GetXaxis()->SetTitle("#chi^2");
 206         histo->GetYaxis()->SetTitle("Entries");
 207         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist));
 208         histo->GetXaxis()->SetTitle("Sector");
 209         histo->GetYaxis()->SetTitle("Number of Tracks");
 210         histo = dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight));
 211         histo->GetXaxis()->SetTitle("Time / 100ns");
 212         histo->GetYaxis()->SetTitle("Average Pulse Height (a. u.)");
 213         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge));
 214         histo->GetXaxis()->SetTitle("Cluster Charge (a.u.)");
 215         histo->GetYaxis()->SetTitle("Entries");
 216         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit));
 217         histo->GetXaxis()->SetTitle("Charge Deposit (a.u.)");
 218         histo->GetYaxis()->SetTitle("Entries");
 219
 220         // Calculate the purity of the trigger clusters
 221         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
 222         TH1F *histoTracks = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
 223         histoTracks->Divide(histo);
 224         Float_t purities[20], val = 0;
 225         TString triggernames[20];
 226         Int_t nTriggerClasses = 0;
 227         for(Int_t ibin = 1; ibin <= histo->GetNbinsX(); ibin++){
 228                 if((val = histoTracks->GetBinContent(ibin))){
 229                         purities[nTriggerClasses] = val;
 230                         triggernames[nTriggerClasses] = histoTracks->GetXaxis()->GetBinLabel(ibin);
 231                         nTriggerClasses++;
 232                 }
 233         }
 234         TH1F *hTriggerInf = new TH1F("fTriggerInf", "Trigger Information", TMath::Max(nTriggerClasses, 1), 0, TMath::Max(nTriggerClasses, 1));
 235         for(Int_t ibin = 1; ibin <= nTriggerClasses; ibin++){
 236                 hTriggerInf->SetBinContent(ibin, purities[ibin-1]);
 237                 hTriggerInf->GetXaxis()->SetBinLabel(ibin, triggernames[ibin-1].Data());
 238         }
 239         hTriggerInf->GetXaxis()->SetTitle("Trigger Cluster");
 240         hTriggerInf->GetYaxis()->SetTitle("Ratio");
 241         hTriggerInf->GetYaxis()->SetRangeUser(0,1);
 242 //      hTriggerInf->SetMarkerColor(kBlue);
 243 //      hTriggerInf->SetMarkerStyle(22);
 244         fContainer->Add(hTriggerInf);
 245         fNRefFigures = 10;
 246         return kTRUE;
 247 }
 248
 249 //_______________________________________________________
 250 void AliTRDcheckDetector::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t *opt){
 251         //
 252         // Setting Reference Figures
 253         //
 254         opt = "pl";
 255         switch(ifig){
 256                 case 0: first = last = kNTracksEventHist;
 257                                                 break;
 258                 case 1: first = last = kNclustersHist;
 259                                                 break;
 260                 case 2: first = last = kNtrackletsHist;
 261                                                 break;
 262                 case 3: first = last = kNclusterTrackletHist;
 263                                                 break;
 264                 case 4: first = last = kChi2;
 265                                                 break;
 266                 case 5: first = last = kNTracksSectorHist;
 267                                                 break;
 268                 case 6: first = last = kPulseHeight;
 269                                                 break;
 270                 case 7: first = last = kClusterCharge;
 271                                                 break;
 272                 case 8: first = last = kChargeDeposit;
 273                                                 break;
 274                 case 9: first = last = kPurity;
 275                                                 opt="bar";
 276                                                 break;
 277                 default: first = last = kNTracksEventHist;
 278                                                 break;
 279         };
 280 }
 281