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