ANALYSIS/AliESDpidCuts.cxx

   1 #include <TCanvas.h>
   2 #include <TClass.h>
   3 #include <TCollection.h>
   4 #include <TDirectory.h>
   5 #include <TH1F.h>
   6 #include <TH1I.h>
   7 #include <TH2I.h>
   8 #include <TMath.h>
   9 #include <TIterator.h>
  10 #include <TString.h>
  11
  12 #include "AliESDtrack.h"
  13 #include "AliESDEvent.h"
  14 #include "AliLog.h"
  15 #include "AliESDpid.h"
  16
  17 #include "AliESDpidCuts.h"
  18
  19 ClassImp(AliESDpidCuts)
  20
  21 const Int_t AliESDpidCuts::kNcuts = 3;
  22
  23 //_____________________________________________________________________
  24 AliESDpidCuts::AliESDpidCuts(const Char_t *name, const Char_t *title):
  25     AliAnalysisCuts(name, title)
  26   , fESDpid(NULL)
  27   , fTPCsigmaCutRequired(0)
  28   , fTOFsigmaCutRequired(0)
  29   , fCutTPCclusterRatio(0.)
  30   , fMinMomentumTOF(0.5)
  31   , fHcutStatistics(NULL)
  32   , fHcutCorrelation(NULL)
  33 {
  34   //
  35   // Default constructor
  36   //
  37
  38   fESDpid = new AliESDpid;
  39   memset(fCutTPCnSigma, 0, sizeof(Float_t) * 2);
  40   memset(fCutTOFnSigma, 0, sizeof(Float_t) * 2);
  41
  42   memset(fHclusterRatio, 0, sizeof(TH1F *) * 2);
  43   memset(fHnSigmaTPC, 0, sizeof(TH1F *) * AliPID::kSPECIES * 2);
  44   memset(fHnSigmaTOF, 0, sizeof(TH1F *) * AliPID::kSPECIES * 2);
  45 }
  46
  47 //_____________________________________________________________________
  48 AliESDpidCuts::AliESDpidCuts(const AliESDpidCuts &ref):
  49     AliAnalysisCuts(ref)
  50   , fESDpid(NULL)
  51   , fTPCsigmaCutRequired(ref.fTPCsigmaCutRequired)
  52   , fTOFsigmaCutRequired(ref.fTOFsigmaCutRequired)
  53   , fCutTPCclusterRatio(ref.fCutTPCclusterRatio)
  54   , fMinMomentumTOF(ref.fMinMomentumTOF)
  55   , fHcutStatistics(NULL)
  56   , fHcutCorrelation(NULL)
  57 {
  58   //
  59   // Copy constructor
  60   //
  61   fESDpid = new AliESDpid(*ref.fESDpid);
  62   memcpy(fCutTPCnSigma, ref.fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
  63   memcpy(fCutTOFnSigma, ref.fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
  64
  65   if(ref.fHcutStatistics) fHcutStatistics = dynamic_cast<TH1I *>(ref.fHcutStatistics->Clone());
  66   if(ref.fHcutCorrelation) fHcutCorrelation = dynamic_cast<TH2I *>(ref.fHcutCorrelation->Clone());
  67   for(Int_t imode = 0; imode < 2; imode++){
  68     if(ref.fHclusterRatio[imode]) fHclusterRatio[imode] = dynamic_cast<TH1F *>(ref.fHclusterRatio[imode]->Clone());
  69     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
  70       if(fHnSigmaTPC[ispec][imode]) fHnSigmaTPC[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());
  71       if(fHnSigmaTOF[ispec][imode]) fHnSigmaTOF[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());
  72     }
  73   }
  74 }
  75
  76 //_____________________________________________________________________
  77 AliESDpidCuts &AliESDpidCuts::operator=(const AliESDpidCuts &ref){
  78   //
  79   // Assignment operator
  80   //
  81   if(this != &ref)
  82     ref.Copy(*this);
  83   return *this;
  84 }
  85
  86 //_____________________________________________________________________
  87 AliESDpidCuts::~AliESDpidCuts(){
  88   //
  89   // Destructor
  90   //
  91   delete fESDpid;
  92
  93   delete fHcutStatistics;
  94   delete fHcutCorrelation;
  95   for(Int_t imode = 0; imode < 2; imode++){
  96     delete fHclusterRatio[imode];
  97     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
  98       delete fHnSigmaTPC[ispec][imode];
  99       delete fHnSigmaTOF[ispec][imode];
 100     }
 101   }
 102 }
 103
 104 //_____________________________________________________________________
 105 Bool_t AliESDpidCuts::IsSelected(TObject *track, TObject *event){
 106   //
 107   // Select Track
 108   //
 109   if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
 110     Char_t errormessage[256];
 111     sprintf(errormessage, "Provided object not an AliESDtrack: Type %s", track->IsA()->GetName());
 112     AliError(errormessage);
 113     return kFALSE;
 114   }
 115   return AcceptTrack(const_cast<const AliESDtrack *>(dynamic_cast<AliESDtrack *>(track)),dynamic_cast<AliESDEvent *>(event));
 116 }
 117
 118 //_____________________________________________________________________
 119 void AliESDpidCuts::Copy(TObject &c) const {
 120   //
 121   // Copy function
 122   //
 123   AliESDpidCuts &target = dynamic_cast<AliESDpidCuts &>(c);
 124
 125   target.fESDpid = new AliESDpid(*fESDpid);
 126
 127   target.fCutTPCclusterRatio = fCutTPCclusterRatio;
 128   target.fMinMomentumTOF = fMinMomentumTOF;
 129
 130   target.fTPCsigmaCutRequired = fTPCsigmaCutRequired;
 131   target.fTOFsigmaCutRequired = fTOFsigmaCutRequired;
 132
 133   if(fHcutStatistics) target.fHcutStatistics = dynamic_cast<TH1I *>(fHcutStatistics->Clone());
 134   if(fHcutCorrelation) target.fHcutCorrelation = dynamic_cast<TH2I *>(fHcutCorrelation->Clone());
 135   for(Int_t imode = 0; imode < 2; imode++){
 136     if(fHclusterRatio[imode]) target.fHclusterRatio[imode] = dynamic_cast<TH1F *>(fHclusterRatio[imode]->Clone());
 137     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 138       if(fHnSigmaTPC[ispec][imode]) target.fHnSigmaTPC[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());
 139       if(fHnSigmaTOF[ispec][imode]) target.fHnSigmaTOF[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTOF[ispec][imode]->Clone());
 140     }
 141   }
 142
 143   memcpy(target.fCutTPCnSigma, fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
 144   memcpy(target.fCutTOFnSigma, fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES);
 145
 146   AliESDpidCuts::Copy(c);
 147 }
 148
 149 //_____________________________________________________________________
 150 Long64_t AliESDpidCuts::Merge(TCollection *coll){
 151   //
 152   // Merge Cut objects
 153   //
 154   if(coll) return 0;
 155   if(coll->IsEmpty()) return 1;
 156   if(!HasHistograms())  return 0;
 157
 158   TIterator *iter = coll->MakeIterator();
 159   TObject *o = NULL;
 160   AliESDpidCuts *ref = NULL;
 161   Int_t counter = 0;
 162   while((o = iter->Next())){
 163     ref = dynamic_cast<AliESDpidCuts *>(o);
 164     if(!ref) continue;
 165     if(!ref->HasHistograms()) continue;
 166
 167     fHcutStatistics->Add(ref->fHcutStatistics);
 168     fHcutCorrelation->Add(ref->fHcutCorrelation);
 169     for(Int_t imode = 0; imode < 2; imode++){
 170       fHclusterRatio[imode]->Add(ref->fHclusterRatio[imode]);
 171       for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 172         fHnSigmaTPC[ispec][imode]->Add(ref->fHnSigmaTPC[ispec][imode]);
 173         fHnSigmaTOF[ispec][imode]->Add(ref->fHnSigmaTOF[ispec][imode]);
 174       }
 175     }
 176     ++counter;
 177   }
 178   return ++counter;
 179 }
 180
 181 //_____________________________________________________________________
 182 void AliESDpidCuts::DefineHistograms(Color_t color){
 183   //
 184   // Swich on QA and create the histograms
 185   //
 186   SetBit(kHasHistograms, kTRUE);
 187   fHcutStatistics = new TH1I("fHcutStatistics", "Cut Statistics", kNcuts, 0, kNcuts);
 188   fHcutStatistics->SetLineColor(color);
 189   fHcutCorrelation = new TH2I("fHcutCorrelation", "Cut Correlation", kNcuts, 0, kNcuts, kNcuts, 0, kNcuts);
 190   TString cutname[kNcuts] = {"TPCclusterRatio", "TPC sigma", "TOF sigma"};
 191   for(Int_t icut = 0; icut < kNcuts; icut++){
 192     fHcutStatistics->GetXaxis()->SetBinLabel(fHcutStatistics->GetXaxis()->GetFirst() + icut, cutname[icut].Data());
 193     fHcutCorrelation->GetXaxis()->SetBinLabel(fHcutCorrelation->GetXaxis()->GetFirst() + icut, cutname[icut].Data());
 194     fHcutCorrelation->GetYaxis()->SetBinLabel(fHcutCorrelation->GetYaxis()->GetFirst() + icut, cutname[icut].Data());
 195   }
 196   Char_t hname[256], htitle[256];
 197   for(Int_t imode = 0; imode < 2; imode++){
 198     sprintf(hname, "fHclusterRatio%s", imode ? "After" : "Before");
 199     sprintf(htitle, "TPC cluster Ratio %s cuts;Ratio;Entries", imode ? "after" : "before");
 200     fHclusterRatio[imode] = new TH1F(hname, htitle, 20, 0., 1.);
 201     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 202       sprintf(hname, "fHnSigma%sTPC%s", AliPID::ParticleName(ispec), imode ? "after" : "before");
 203       sprintf(htitle, "TPC sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");
 204       fHnSigmaTPC[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);
 205       sprintf(hname, "fHnSigma%sTOF%s", AliPID::ParticleName(ispec), imode ? "after" : "before");
 206       sprintf(htitle, "TOF sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");
 207       fHnSigmaTOF[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);
 208     }
 209   }
 210 }
 211
 212 //_____________________________________________________________________
 213 Bool_t AliESDpidCuts::AcceptTrack(const AliESDtrack *track, const AliESDEvent *event){
 214   //
 215   // Check whether the tracks survived the cuts
 216   //
 217   enum{
 218     kCutClusterRatioTPC,
 219     kCutNsigmaTPC,
 220     kCutNsigmaTOF
 221   };
 222   Long64_t cutRequired=0, cutFullfiled = 0;
 223   if(fTOFsigmaCutRequired && event == 0)
 224     AliError("No event pointer. Need event pointer for T0 for TOF cut");
 225   Double_t clusterRatio = track->GetTPCNclsF() ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;
 226   if(fCutTPCclusterRatio > 0.){
 227     SETBIT(cutRequired, kCutClusterRatioTPC);
 228     if(clusterRatio >= fCutTPCclusterRatio)
 229       SETBIT(cutFullfiled, kCutClusterRatioTPC);
 230   }
 231   // check TPC nSigma cut
 232   Float_t nsigmaTPC[AliPID::kSPECIES], nsigma;   // need all sigmas for QA plotting
 233   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 234     nsigmaTPC[ispec] = nsigma = fESDpid->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);
 235     if(!(fTPCsigmaCutRequired & 1 << ispec)) continue;
 236     SETBIT(cutRequired, kCutNsigmaTPC); // We found at least one species where the n-Sigma Cut is required
 237     if(nsigma >= fCutTPCnSigma[2*ispec] && nsigma <= fCutTPCnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTPC);    // Fullfiled for at least one species
 238   }
 239   // check TOF nSigma cut
 240   Float_t nsigmaTOF[AliPID::kSPECIES];    // see above
 241   Bool_t hasTOFpid = track->GetStatus() & AliESDtrack::kTOFpid; // only apply TOF n-sigma cut when PID Status Bit is set
 242   Double_t times[AliPID::kSPECIES];
 243   track->GetIntegratedTimes(times);
 244   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 245
 246     if(hasTOFpid) nsigmaTOF[ispec] = nsigma = fESDpid->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec, event->GetT0());
 247     if(!(fTOFsigmaCutRequired && 1 << ispec)) continue;
 248     SETBIT(cutRequired, kCutNsigmaTOF);
 249     if(track->GetOuterParam()->P() >= fMinMomentumTOF){
 250       if(hasTOFpid && nsigma <= fCutTOFnSigma[2*ispec] && nsigma >= fCutTOFnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTOF);
 251     }
 252   }
 253
 254   // Fill Histograms before cuts
 255   if(HasHistograms()){
 256     fHclusterRatio[0]->Fill(clusterRatio);
 257     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 258       fHnSigmaTPC[ispec][0]->Fill(nsigmaTPC[ispec]);
 259       if(hasTOFpid) fHnSigmaTOF[ispec][0]->Fill(nsigmaTOF[ispec]);
 260     }
 261   }
 262   if(cutRequired != cutFullfiled){
 263     // Fill cut statistics
 264     if(HasHistograms()){
 265       for(Int_t icut = 0; icut < kNcuts; icut++){
 266         if(TESTBIT(cutRequired, icut) && !TESTBIT(cutFullfiled, icut)){
 267           // cut not fullfiled
 268           fHcutStatistics->Fill(icut);
 269           for(Int_t jcut = 0; jcut <= icut; jcut++)
 270             if(TESTBIT(cutRequired, jcut) && !TESTBIT(cutFullfiled, jcut)) fHcutCorrelation->Fill(jcut, icut);
 271         }
 272       }
 273     }
 274     return kFALSE;    // At least one cut is not fullfiled
 275   }
 276
 277   // Fill Histograms after cuts
 278   if(HasHistograms()){
 279     fHclusterRatio[1]->Fill(clusterRatio);
 280     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 281       fHnSigmaTPC[ispec][1]->Fill(nsigmaTPC[ispec]);
 282       if(hasTOFpid) fHnSigmaTOF[ispec][1]->Fill(nsigmaTOF[ispec]);
 283     }
 284   }
 285
 286   return kTRUE;
 287 }
 288
 289 //_____________________________________________________________________
 290 void AliESDpidCuts::SaveHistograms(const Char_t * location){
 291   //
 292   // Save the histograms to a file
 293   //
 294   if(!HasHistograms()){
 295     AliError("Histograms not on - Exiting");
 296     return;
 297   }
 298   if(!location) location = GetName();
 299   gDirectory->mkdir(location);
 300   gDirectory->cd(location);
 301   fHcutStatistics->Write();
 302   fHcutCorrelation->Write();
 303
 304   gDirectory->mkdir("before_cuts");
 305   gDirectory->mkdir("after_cuts");
 306
 307   gDirectory->cd("before_cuts");
 308   fHclusterRatio[0]->Write();
 309   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 310     fHnSigmaTPC[ispec][0]->Write();
 311     fHnSigmaTOF[ispec][0]->Write();
 312   }
 313
 314   gDirectory->cd("../after_cuts");
 315   fHclusterRatio[1]->Write();
 316   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 317     fHnSigmaTPC[ispec][1]->Write();
 318     fHnSigmaTOF[ispec][1]->Write();
 319   }
 320
 321   gDirectory->cd("..");
 322 }
 323
 324 //_____________________________________________________________________
 325 void AliESDpidCuts::DrawHistograms(){
 326   //
 327   // Draw the Histograms
 328   //
 329   TCanvas *stat = new TCanvas("cutStat", "Cut Statistics", 640, 480);
 330   stat->cd();
 331   fHcutStatistics->SetStats(kFALSE);
 332   fHcutStatistics->Draw();
 333   stat->SaveAs(Form("%s_%s.gif", GetName(), stat->GetName()));
 334
 335   TCanvas *correl = new TCanvas("cutCorrelation", "Cut Correlation", 640, 480);
 336   correl->cd();
 337   fHcutCorrelation->SetStats(kFALSE);
 338   fHcutCorrelation->Draw("colz");
 339   correl->SaveAs(Form("%s_%s.gif", GetName(), correl->GetName()));
 340
 341   TCanvas *cRatio = new TCanvas("ClusterRatioTPC", "TPC cluster Ratio", 640, 480);
 342   cRatio->cd();
 343   fHclusterRatio[0]->SetLineColor(kRed);
 344   fHclusterRatio[0]->SetStats(kFALSE);
 345   fHclusterRatio[0]->Draw();
 346   fHclusterRatio[1]->SetLineColor(kBlue);
 347   fHclusterRatio[1]->SetStats(kFALSE);
 348   fHclusterRatio[1]->Draw("same");
 349   cRatio->SaveAs(Form("%s_%s.gif",  GetName(), cRatio->GetName()));
 350
 351   TCanvas *cNsigTPC = new TCanvas("NsigmaTPC", "TPC n-sigma", 640, 480);
 352   cNsigTPC->Divide(3,2);
 353   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 354     cNsigTPC->cd(ispec + 1);
 355     fHnSigmaTPC[ispec][0]->SetLineColor(kRed);
 356     fHnSigmaTPC[ispec][0]->SetStats(kFALSE);
 357     fHnSigmaTPC[ispec][0]->Draw();
 358     fHnSigmaTPC[ispec][1]->SetLineColor(kBlue);
 359     fHnSigmaTPC[ispec][1]->SetStats(kFALSE);
 360     fHnSigmaTPC[ispec][1]->Draw("same");
 361   }
 362   cNsigTPC->SaveAs(Form("%s_%s.gif", GetName(), cNsigTPC->GetName()));
 363
 364   TCanvas *cNsigTOF = new TCanvas("NsigmaTOF", "TOF n-sigma", 640, 480);
 365   cNsigTOF->Divide(3,2);
 366   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 367     cNsigTOF->cd(ispec + 1);
 368     fHnSigmaTOF[ispec][0]->SetLineColor(kRed);
 369     fHnSigmaTOF[ispec][0]->SetStats(kFALSE);
 370     fHnSigmaTOF[ispec][0]->Draw();
 371     fHnSigmaTOF[ispec][1]->SetLineColor(kBlue);
 372     fHnSigmaTOF[ispec][1]->SetStats(kFALSE);
 373     fHnSigmaTOF[ispec][1]->Draw("same");
 374   }
 375   cNsigTOF->SaveAs(Form("%s_%s.gif", GetName(), cNsigTOF->GetName()));
 376 }
 377