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     Init();\r
 251     if (!fPIDresponse)\r
 252       AliError("PID Response not available");\r
 253     return 0;\r
 254   }\r
 255   enum{\r
 256     kCutClusterRatioTPC,\r
 257     kCutNsigmaTPC,\r
 258     kCutNsigmaTOF\r
 259   };\r
 260   Long64_t cutRequired=0, cutFullfiled = 0;\r
 261   if(fTOFsigmaCutRequired && event == 0)  {\r
 262     AliError("No event pointer. Need event pointer for T0 for TOF cut");\r
 263     return (0);\r
 264   }\r
 265   Double_t clusterRatio = track->GetTPCNclsF() ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;\r
 266   if(fCutTPCclusterRatio > 0.){\r
 267     SETBIT(cutRequired, kCutClusterRatioTPC);\r
 268     if(clusterRatio >= fCutTPCclusterRatio) \r
 269       SETBIT(cutFullfiled, kCutClusterRatioTPC);\r
 270   }\r
 271   // check TPC nSigma cut\r
 272   Float_t nsigmaTPC[AliPID::kSPECIES], nsigma;   // need all sigmas for QA plotting\r
 273   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
 274     nsigmaTPC[ispec] = nsigma = fPIDresponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
 275     if(!(fTPCsigmaCutRequired & 1 << ispec)) continue;\r
 276     SETBIT(cutRequired, kCutNsigmaTPC); // We found at least one species where the n-Sigma Cut is required\r
 277     if(nsigma >= fCutTPCnSigma[2*ispec] && nsigma <= fCutTPCnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTPC);    // Fullfiled for at least one species\r
 278   }\r
 279   // check TOF nSigma cut\r
 280   Float_t nsigmaTOF[AliPID::kSPECIES];    // see above\r
 281   Bool_t hasTOFpid = track->GetStatus() & AliESDtrack::kTOFpid; // only apply TOF n-sigma cut when PID Status Bit is set\r
 282   Double_t times[AliPID::kSPECIES];\r
 283   track->GetIntegratedTimes(times);\r
 284   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
 285     \r
 286     if(hasTOFpid && event) nsigmaTOF[ispec] = nsigma = fPIDresponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);\r
 287     if(!(fTOFsigmaCutRequired && 1 << ispec)) continue;\r
 288     SETBIT(cutRequired, kCutNsigmaTOF);\r
 289     if(track->GetOuterParam()->P() >= fMinMomentumTOF){\r
 290       if(hasTOFpid && nsigma <= fCutTOFnSigma[2*ispec] && nsigma >= fCutTOFnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTOF);\r
 291     }\r
 292   }\r
 293 \r
 294   // Fill Histograms before cuts\r
 295   if(HasHistograms()){\r
 296     fHclusterRatio[0]->Fill(clusterRatio);\r
 297     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
 298       fHnSigmaTPC[ispec][0]->Fill(nsigmaTPC[ispec]);\r
 299       if(hasTOFpid) fHnSigmaTOF[ispec][0]->Fill(nsigmaTOF[ispec]);\r
 300     }\r
 301   }\r
 302   if(cutRequired != cutFullfiled){\r
 303     // Fill cut statistics\r
 304     if(HasHistograms()){\r
 305       for(Int_t icut = 0; icut < kNcuts; icut++){\r
 306         if(TESTBIT(cutRequired, icut) && !TESTBIT(cutFullfiled, icut)){\r
 307           // cut not fullfiled\r
 308           fHcutStatistics->Fill(icut);\r
 309           for(Int_t jcut = 0; jcut <= icut; jcut++)\r
 310             if(TESTBIT(cutRequired, jcut) && !TESTBIT(cutFullfiled, jcut)) fHcutCorrelation->Fill(jcut, icut);\r
 311         }\r
 312       }\r
 313     }\r
 314     return kFALSE;    // At least one cut is not fullfiled\r
 315   }\r
 316 \r
 317   // Fill Histograms after cuts\r
 318   if(HasHistograms()){\r
 319     fHclusterRatio[1]->Fill(clusterRatio);\r
 320     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
 321       fHnSigmaTPC[ispec][1]->Fill(nsigmaTPC[ispec]);\r
 322       if(hasTOFpid) fHnSigmaTOF[ispec][1]->Fill(nsigmaTOF[ispec]);\r
 323     }\r
 324   }\r
 325 \r
 326   return kTRUE;\r
 327 }\r
 328 \r
 329 //_____________________________________________________________________\r
 330 void AliESDpidCuts::SaveHistograms(const Char_t * location){\r
 331   //\r
 332   // Save the histograms to a file\r
 333   //\r
 334   if(!HasHistograms()){\r
 335     AliError("Histograms not on - Exiting");\r
 336     return;\r
 337   }\r
 338   if(!location) location = GetName();\r
 339   gDirectory->mkdir(location);\r
 340   gDirectory->cd(location);\r
 341   fHcutStatistics->Write();\r
 342   fHcutCorrelation->Write();\r
 343 \r
 344   gDirectory->mkdir("before_cuts");\r
 345   gDirectory->mkdir("after_cuts");\r
 346 \r
 347   gDirectory->cd("before_cuts");\r
 348   fHclusterRatio[0]->Write();\r
 349   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
 350     fHnSigmaTPC[ispec][0]->Write();\r
 351     fHnSigmaTOF[ispec][0]->Write();\r
 352   }\r
 353 \r
 354   gDirectory->cd("../after_cuts");\r
 355   fHclusterRatio[1]->Write();\r
 356   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
 357     fHnSigmaTPC[ispec][1]->Write();\r
 358     fHnSigmaTOF[ispec][1]->Write();\r
 359   }\r
 360 \r
 361   gDirectory->cd("..");\r
 362 }\r
 363 \r
 364 //_____________________________________________________________________\r
 365 void AliESDpidCuts::DrawHistograms(){\r
 366   //\r
 367   // Draw the Histograms\r
 368   //\r
 369   TCanvas *stat = new TCanvas("cutStat", "Cut Statistics", 640, 480);\r
 370   stat->cd();\r
 371   fHcutStatistics->SetStats(kFALSE);\r
 372   fHcutStatistics->Draw();\r
 373   stat->SaveAs(Form("%s_%s.gif", GetName(), stat->GetName()));\r
 374 \r
 375   TCanvas *correl = new TCanvas("cutCorrelation", "Cut Correlation", 640, 480);\r
 376   correl->cd();\r
 377   fHcutCorrelation->SetStats(kFALSE);\r
 378   fHcutCorrelation->Draw("colz");\r
 379   correl->SaveAs(Form("%s_%s.gif", GetName(), correl->GetName()));\r
 380 \r
 381   TCanvas *cRatio = new TCanvas("ClusterRatioTPC", "TPC cluster Ratio", 640, 480);\r
 382   cRatio->cd();\r
 383   fHclusterRatio[0]->SetLineColor(kRed);\r
 384   fHclusterRatio[0]->SetStats(kFALSE);\r
 385   fHclusterRatio[0]->Draw();\r
 386   fHclusterRatio[1]->SetLineColor(kBlue);\r
 387   fHclusterRatio[1]->SetStats(kFALSE);\r
 388   fHclusterRatio[1]->Draw("same");\r
 389   cRatio->SaveAs(Form("%s_%s.gif",  GetName(), cRatio->GetName()));\r
 390 \r
 391   TCanvas *cNsigTPC = new TCanvas("NsigmaTPC", "TPC n-sigma", 640, 480);\r
 392   cNsigTPC->Divide(3,2);\r
 393   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
 394     cNsigTPC->cd(ispec + 1);\r
 395     fHnSigmaTPC[ispec][0]->SetLineColor(kRed);\r
 396     fHnSigmaTPC[ispec][0]->SetStats(kFALSE);\r
 397     fHnSigmaTPC[ispec][0]->Draw();\r
 398     fHnSigmaTPC[ispec][1]->SetLineColor(kBlue);\r
 399     fHnSigmaTPC[ispec][1]->SetStats(kFALSE);\r
 400     fHnSigmaTPC[ispec][1]->Draw("same");\r
 401   }\r
 402   cNsigTPC->SaveAs(Form("%s_%s.gif", GetName(), cNsigTPC->GetName()));\r
 403 \r
 404   TCanvas *cNsigTOF = new TCanvas("NsigmaTOF", "TOF n-sigma", 640, 480);\r
 405   cNsigTOF->Divide(3,2);\r
 406   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
 407     cNsigTOF->cd(ispec + 1);\r
 408     fHnSigmaTOF[ispec][0]->SetLineColor(kRed);\r
 409     fHnSigmaTOF[ispec][0]->SetStats(kFALSE);\r
 410     fHnSigmaTOF[ispec][0]->Draw();\r
 411     fHnSigmaTOF[ispec][1]->SetLineColor(kBlue);\r
 412     fHnSigmaTOF[ispec][1]->SetStats(kFALSE);\r
 413     fHnSigmaTOF[ispec][1]->Draw("same");\r
 414   }\r
 415   cNsigTOF->SaveAs(Form("%s_%s.gif", GetName(), cNsigTOF->GetName()));\r
 416 }\r
 417 \r