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