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           Int_t crossing = fTracklet->GetNChange();
 151           if(fDebugLevel > 3){
 152             (*fDebugStream) << "ChargeDeposit"
 153               << "Detector="  << detector
 154               << "Sector="    << sector
 155               << "nclusters=" << nClustersTracklet
 156               << "crossing="  << crossing
 157               << "QT="        << Qtot
 158               << "\n";
 159           }
 160         }
 161       }
 162       dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist))->Fill(sector);
 163     }
 164     nTracks++;
 165   }
 166   if(nTracks){
 167     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks))->Fill(triggermask);
 168     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist))->Fill(nTracks);
 169   }
 170   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
 171     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
 172     // also set the label for both histograms
 173     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
 174     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 175     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
 176     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 177   }
 178   PostData(0, fContainer);
 179 }
 180
 181 //_______________________________________________________
 182 void AliTRDcheckDetector::Terminate(Option_t *){
 183   //
 184   // Terminate function
 185   //
 186 }
 187
 188 //_______________________________________________________
 189 Bool_t AliTRDcheckDetector::PostProcess(){
 190   //
 191   // Do Postprocessing (for the moment set the number of Reference histograms)
 192   //
 193
 194   TH1 * histo = 0x0;
 195   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist));
 196   histo->GetXaxis()->SetTitle("Number of Tracks");
 197   histo->GetYaxis()->SetTitle("Events");
 198   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist));
 199   histo->GetXaxis()->SetTitle("Number of Clusters");
 200   histo->GetYaxis()->SetTitle("Entries");
 201   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist));
 202   histo->GetXaxis()->SetTitle("Number of Tracklets");
 203   histo->GetYaxis()->SetTitle("Entries");
 204   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist));
 205   histo->GetXaxis()->SetTitle("Number of Clusters");
 206   histo->GetYaxis()->SetTitle("Entries");
 207   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2));
 208   histo->GetXaxis()->SetTitle("#chi^2");
 209   histo->GetYaxis()->SetTitle("Entries");
 210   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist));
 211   histo->GetXaxis()->SetTitle("Sector");
 212   histo->GetYaxis()->SetTitle("Number of Tracks");
 213   histo = dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight));
 214   histo->GetXaxis()->SetTitle("Time / 100ns");
 215   histo->GetYaxis()->SetTitle("Average Pulse Height (a. u.)");
 216   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge));
 217   histo->GetXaxis()->SetTitle("Cluster Charge (a.u.)");
 218   histo->GetYaxis()->SetTitle("Entries");
 219   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit));
 220   histo->GetXaxis()->SetTitle("Charge Deposit (a.u.)");
 221   histo->GetYaxis()->SetTitle("Entries");
 222
 223   // Calculate the purity of the trigger clusters
 224   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
 225   TH1F *histoTracks = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
 226   histoTracks->Divide(histo);
 227   Float_t purities[20], val = 0;
 228   TString triggernames[20];
 229   Int_t nTriggerClasses = 0;
 230   for(Int_t ibin = 1; ibin <= histo->GetNbinsX(); ibin++){
 231     if((val = histoTracks->GetBinContent(ibin))){
 232       purities[nTriggerClasses] = val;
 233       triggernames[nTriggerClasses] = histoTracks->GetXaxis()->GetBinLabel(ibin);
 234       nTriggerClasses++;
 235     }
 236   }
 237   TH1F *hTriggerInf = new TH1F("fTriggerInf", "Trigger Information", TMath::Max(nTriggerClasses, 1), 0, TMath::Max(nTriggerClasses, 1));
 238   for(Int_t ibin = 1; ibin <= nTriggerClasses; ibin++){
 239     hTriggerInf->SetBinContent(ibin, purities[ibin-1]);
 240     hTriggerInf->GetXaxis()->SetBinLabel(ibin, triggernames[ibin-1].Data());
 241   }
 242   hTriggerInf->GetXaxis()->SetTitle("Trigger Cluster");
 243   hTriggerInf->GetYaxis()->SetTitle("Ratio");
 244   hTriggerInf->GetYaxis()->SetRangeUser(0,1);
 245 //      hTriggerInf->SetMarkerColor(kBlue);
 246 //      hTriggerInf->SetMarkerStyle(22);
 247   fContainer->Add(hTriggerInf);
 248   fNRefFigures = 10;
 249   return kTRUE;
 250 }
 251
 252 //_______________________________________________________
 253 void AliTRDcheckDetector::GetRefFigure(Int_t ifig){
 254   //
 255   // Setting Reference Figures
 256   //
 257   TH1 *h = 0x0;
 258   switch(ifig){
 259   case 0:
 260     ((TH1F*)fContainer->At(kNTracksEventHist))->Draw("pl");
 261     break;
 262   case 1:
 263     ((TH1F*)fContainer->At(kNclustersHist))->Draw("pl");
 264     break;
 265   case 2:
 266     h = (TH1F*)fContainer->At(kNtrackletsHist);
 267     if(!h->GetEntries()) break;
 268     h->Scale(100./h->Integral());
 269     h->GetXaxis()->SetRangeUser(.5, 6.5);
 270     h->SetFillColor(kGreen);
 271     h->SetBarOffset(.2);
 272     h->SetBarWidth(.6);
 273     h->Draw("bar1");
 274     break;
 275   case 3:
 276     ((TH1F*)fContainer->At(kNclusterTrackletHist))->Draw("pc");
 277     break;
 278   case 4:
 279     ((TH1F*)fContainer->At(kChi2))->Draw("");
 280     break;
 281   case 5:
 282     h = (TH1F*)fContainer->At(kNTracksSectorHist);
 283     if(!h->GetEntries()) break;
 284     h->Scale(100./h->Integral());
 285     h->SetFillColor(kGreen);
 286     h->SetBarOffset(.2);
 287     h->SetBarWidth(.6);
 288     h->Draw("bar1");
 289     break;
 290   case 6:
 291     h = (TH1F*)fContainer->At(kPulseHeight);
 292     h->SetMarkerStyle(24);
 293     h->Draw("e1");
 294     break;
 295   case 7:
 296     ((TH1F*)fContainer->At(kClusterCharge))->Draw("c");
 297     break;
 298   case 8:
 299     ((TH1F*)fContainer->At(kChargeDeposit))->Draw("c");
 300     break;
 301   case 9:
 302     h=(TH1F*)fContainer->At(kPurity);
 303     h->SetBarOffset(.2);
 304     h->SetBarWidth(.6);
 305     h->SetFillColor(kGreen);
 306     h->Draw("bar1");
 307     break;
 308   default:
 309     ((TH1F*)fContainer->At(kNTracksEventHist))->Draw("pl");
 310     break;
 311   }
 312 }
 313