TRD/qaRec/AliTRDcheckDetector.cxx

   1 #include <TFile.h>
   2 #include <TH1F.h>
   3 #include <TMath.h>
   4 #include <TObjArray.h>
   5 #include <TProfile.h>
   6 #include <TROOT.h>
   7
   8 #include "AliTRDcluster.h"
   9 #include "AliTRDseedV1.h"
  10 #include "AliTRDtrackV1.h"
  11 #include "TTreeStream.h"
  12
  13 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  14 #include "AliTRDcheckDetector.h"
  15
  16 #include <cstdio>
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  Reconstruction QA                                                     //
  21 //                                                                        //
  22 //  Task doing basic checks for tracking and detector performance         //
  23 //                                                                        //
  24 //  Authors:                                                              //
  25 //    Anton Andronic <A.Andronic@gsi.de>                                  //
  26 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  27 //                                                                        //
  28 ////////////////////////////////////////////////////////////////////////////
  29
  30 //_______________________________________________________
  31 AliTRDcheckDetector::AliTRDcheckDetector():
  32   AliTRDrecoTask("DetChecker", "Basic Detector Checker")
  33   ,fPHSdetector(0x0)
  34   ,fPHSsector(0x0)
  35   ,fQCLdetector(0x0)
  36   ,fQCLsector(0x0)
  37   ,fQTdetector(0x0)
  38   ,fQTsector(0x0)
  39 {
  40   //
  41   // Default constructor
  42   //
  43 }
  44
  45 //_______________________________________________________
  46 AliTRDcheckDetector::~AliTRDcheckDetector(){
  47   //
  48   // Destructor
  49   //
  50   if(fPHSdetector) delete fPHSdetector;
  51   if(fPHSsector) delete fPHSsector;
  52   if(fQCLdetector) delete fQCLdetector;
  53   if(fQCLsector) delete fQCLsector;
  54   if(fQTdetector) delete fQTdetector;
  55   if(fQTsector) delete fQTsector;
  56 }
  57
  58 //_______________________________________________________
  59 void AliTRDcheckDetector::CreateOutputObjects(){
  60   //
  61   // Create Output Objects
  62   //
  63   fContainer = new TObjArray();
  64   TH1F *hNtrks = new TH1F("hNtrks", "Number of Tracks per event", 100, 0, 100);
  65   TH1F *hNcls = new TH1F("hNcls", "Nr. of clusters per track", 30, 0, 180);
  66   TH1F *hNtls = new TH1F("hNtls", "Nr. tracklets per track", 10, 0, 10);
  67   TH1F *hNclTls = new TH1F("hNclTls","Mean Number of clusters per tracklet", 30, 0, 30);
  68   TH1F *hChi2 = new TH1F("hChi2", "Chi2", 200, 0, 20);
  69   TH1F *hChi2n = new TH1F("hChi2n", "Norm. Chi2 (tracklets)", 50, 0, 5);
  70   TH1F *hSM = new TH1F("hSM", "Track Counts in Supermodule", 18, 0., 18.);
  71   TH1F *hQcl = new TH1F("hQcl1", "Cluster charge", 200, -1200, 1200);
  72   TProfile *hPH = new TProfile("hPH", "Average PH", 31, -0.5, 30.5);
  73   TH1F * hQT = new TH1F("hQT", "Total Charge Deposit", 6000, 0, 6000);
  74   // Register Histograms
  75   fContainer->Add(hNtrks);
  76   fContainer->Add(hNcls);
  77   fContainer->Add(hNtls);
  78   fContainer->Add(hNclTls);
  79   fContainer->Add(hChi2);
  80   fContainer->Add(hChi2n);
  81   fContainer->Add(hSM);
  82   fContainer->Add(hPH);
  83   fContainer->Add(hQcl);
  84   fContainer->Add(hQT);
  85
  86   // Initialise the PHS, cluster charge and total charge deposit histograms for each detector
  87   fPHSdetector = new TObjArray();
  88   fPHSdetector->SetName("fPHSdetector");
  89   fQCLdetector = new TObjArray();
  90   fQCLdetector->SetName("fQCLdetector");
  91   fQTdetector = new TObjArray();
  92   fQTdetector->SetName("fQTdetector");
  93   for(Int_t idet = 0; idet < kNDetectors; idet++){
  94     fPHSdetector->Add(new TProfile(Form("hPHSdet%d", idet), Form("Average Pulse Height in Detector %d", idet), 31, -0.5, 30.5));
  95     fQCLdetector->Add(new TH1F(Form("hQd%d",idet),Form("Qd%d",idet), 100,0,5000));
  96     fQTdetector->Add(new TH1F(Form("hQTd%d", idet), Form("Qtotal%d", idet), 6000, 0, 6000));
  97   }
  98   fContainer->Add(fPHSdetector);
  99   fContainer->Add(fQCLdetector);
 100   fContainer->Add(fQTdetector);
 101
 102   // Initialise the PHS, cluster charge and total charge deposit histograms for each sector
 103   fPHSsector = new TObjArray();
 104   fPHSsector->SetName("fPHSsector");
 105   fQCLsector = new TObjArray();
 106   fQCLsector->SetName("fQCLsector");
 107   fQTsector = new TObjArray();
 108   fQTsector->SetName("fQTsector");
 109   for(Int_t isec = 0; isec < kNSectors; isec++){
 110     fPHSsector->Add(new TProfile(Form("hPHSsec%d", isec), Form("Average Pulse Height in Sector %d", isec), 31, -0.5, 30.5));
 111     fQCLsector->Add(new TH1F(Form("hQd%d",isec),Form("Qd%d",isec), 100,0,5000));
 112     fQTsector->Add(new TH1F(Form("hQTd%d", isec), Form("Qtotal%d", isec), 6000, 0, 6000));
 113   }
 114   fContainer->Add(fPHSsector);
 115   fContainer->Add(fQCLsector);
 116   fContainer->Add(fQTsector);
 117 }
 118
 119 //_______________________________________________________
 120 void AliTRDcheckDetector::Exec(Option_t *){
 121   //
 122   // Execution function
 123   // Filling TRD quality histos
 124   //
 125   Int_t nTracks = 0;            // Count the number of tracks per event
 126   AliTRDtrackInfo *fTrackInfo = 0x0;
 127   AliTRDtrackV1 *fTRDtrack = 0x0;
 128   AliTRDseedV1 *fTracklet = 0x0;
 129   AliTRDcluster *fTRDcluster = 0x0;
 130   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
 131     fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
 132     if(!fTrackInfo || !(fTRDtrack = fTrackInfo->GetTRDtrack())) continue;
 133     Int_t nclusters = fTRDtrack->GetNumberOfClusters();
 134     Int_t ntracklets = fTRDtrack->GetNumberOfTracklets();
 135     Float_t chi2 = fTRDtrack->GetChi2();
 136     // Fill Histograms
 137     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist))->Fill(nclusters);
 138     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist))->Fill(ntracklets);
 139     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2))->Fill(chi2);
 140     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2Normalized))->Fill(chi2/static_cast<Float_t>(ntracklets));
 141     // now loop over single tracklets
 142     if(ntracklets > 2){
 143       Int_t sector = -1;
 144       for(Int_t ilayer = 0; ilayer < kNLayers; ilayer++){
 145         if(!(fTracklet = fTRDtrack->GetTracklet(ilayer)) || !fTracklet->IsOK()) continue;
 146         Float_t nClustersTracklet = fTracklet->GetN();
 147         dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist))->Fill(nClustersTracklet);
 148         if(nClustersTracklet){
 149           Float_t Qtot = 0;
 150           Int_t detector = -1;
 151           for(Int_t itb = 0; itb < kNTimebins; itb++){
 152             if(!(fTRDcluster = fTracklet->GetClusters(itb))) continue;
 153             Int_t localtime        = fTRDcluster->GetLocalTimeBin();
 154             Double_t clusterCharge = fTRDcluster->GetQ();
 155             detector               = fTRDcluster->GetDetector();
 156             sector                 = static_cast<Int_t>(detector/kNDetectorsSector);
 157             Double_t absolute_charge = TMath::Abs(clusterCharge);
 158             Qtot += absolute_charge;
 159             // Fill Histograms
 160             // 1st. PHS
 161             dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight))->Fill(localtime, absolute_charge);
 162             dynamic_cast<TProfile *>(fPHSdetector->UncheckedAt(detector))->Fill(localtime, absolute_charge);
 163             dynamic_cast<TProfile *>(fPHSsector->UncheckedAt(sector))->Fill(localtime, absolute_charge);
 164             // 2nd. cluster charge
 165             dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge))->Fill(absolute_charge);
 166             dynamic_cast<TH1F *>(fQCLdetector->UncheckedAt(detector))->Fill(absolute_charge);
 167             dynamic_cast<TH1F *>(fQCLsector->UncheckedAt(sector))->Fill(absolute_charge);
 168           }
 169           // Fill total cluster charge histograms
 170           dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit))->Fill(Qtot);
 171           dynamic_cast<TH1F *>(fQTdetector->UncheckedAt(detector))->Fill(Qtot);
 172           dynamic_cast<TH1F *>(fQTsector->UncheckedAt(sector))->Fill(Qtot);
 173         }
 174       }
 175       dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist))->Fill(sector);
 176     }
 177     nTracks++;
 178   }
 179   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist))->Fill(nTracks);
 180   PostData(0, fContainer);
 181 }
 182
 183 //_______________________________________________________
 184 void AliTRDcheckDetector::Terminate(Option_t *){
 185   //
 186   // Terminate function
 187   //
 188 }
 189
 190 //_______________________________________________________
 191 Bool_t AliTRDcheckDetector::PostProcess(){
 192         //
 193         // Do Postprocessing (for the moment set the number of Reference histograms)
 194         //
 195
 196         TH1 * histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist));
 197         histo->GetXaxis()->SetTitle("Number of Tracks");
 198         histo->GetYaxis()->SetTitle("Events");
 199         // Calculate percentage of events containing tracks
 200         Int_t entries = (Int_t)histo->GetEntries();
 201   Float_t percTracks = entries > 0 ? (entries - histo->GetBinContent(1))/Float_t(entries) : 0.;
 202         printf("Percentage of Events containing TRD Tracks: %6.3f %%\n", 1.E2*percTracks);
 203         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist));
 204         histo->GetXaxis()->SetTitle("Number of Clusters");
 205         histo->GetYaxis()->SetTitle("Entries");
 206         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist));
 207         histo->GetXaxis()->SetTitle("Number of Tracklets");
 208         histo->GetYaxis()->SetTitle("Entries");
 209         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist));
 210         histo->GetXaxis()->SetTitle("Number of Clusters");
 211         histo->GetYaxis()->SetTitle("Entries");
 212         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2));
 213         histo->GetXaxis()->SetTitle("#chi^2");
 214         histo->GetYaxis()->SetTitle("Entries");
 215         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist));
 216         histo->GetXaxis()->SetTitle("Sector");
 217         histo->GetYaxis()->SetTitle("Number of Tracks");
 218         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist));
 219         histo->GetXaxis()->SetTitle("Sector");
 220         histo->GetYaxis()->SetTitle("Number of Tracks");
 221         histo = dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight));
 222         histo->GetXaxis()->SetTitle("Time / 100ns");
 223         histo->GetYaxis()->SetTitle("Average Pulse Height (a. u.)");
 224         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge));
 225         histo->GetXaxis()->SetTitle("Cluster Charge (a.u.)");
 226         histo->GetYaxis()->SetTitle("Entries");
 227         histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit));
 228         histo->GetXaxis()->SetTitle("Charge Deposit (a.u.)");
 229         histo->GetYaxis()->SetTitle("Entries");
 230         fNRefFigures = 9;
 231         return kTRUE;
 232 }
 233
 234 //_______________________________________________________
 235 void AliTRDcheckDetector::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last){
 236         //
 237         // Setting Reference Figures
 238         //
 239         switch(ifig){
 240                 case 0: first = last = 0;
 241                                                 break;
 242                 case 1: first = last = 1;
 243                                                 break;
 244                 case 2: first = last = 2;
 245                                                 break;
 246                 case 3: first = last = 3;
 247                                                 break;
 248                 case 4: first = last = 4;
 249                                                 break;
 250                 case 5: first = last = 6;
 251                                                 break;
 252                 case 6: first = last = 7;
 253                                                 break;
 254                 case 7: first = last = 8;
 255                                                 break;
 256                 case 8: first = last = 9;
 257                                                 break;
 258                 default: first = last = 0;
 259                                                 break;
 260         };
 261 }
 262