ANALYSIS/AliESDpidCuts.cxx

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