ANALYSIS/AliESDpidCuts.cxx

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