PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraBothTrackCuts.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 //-----------------------------------------------------------------
  18 //         AliSpectraBothTrackCuts class
  19 //-----------------------------------------------------------------
  20
  21 #include "TChain.h"
  22 #include "TTree.h"
  23 #include "TLegend.h"
  24 #include "TH1F.h"
  25 #include "TH1I.h"
  26 #include "TH2F.h"
  27 #include "TCanvas.h"
  28 #include "AliAnalysisTask.h"
  29 #include "AliAnalysisManager.h"
  30 #include "AliAODTrack.h"
  31 #include "AliVTrack.h"
  32 #include "AliExternalTrackParam.h"
  33 #include "AliAODMCParticle.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAODInputHandler.h"
  36 #include "AliAnalysisTaskESDfilter.h"
  37 #include "AliAnalysisDataContainer.h"
  38 #include "AliSpectraBothTrackCuts.h"
  39 //#include "AliSpectraBothHistoManager.h"
  40 #include <iostream>
  41
  42 using namespace std;
  43
  44 const char * AliSpectraBothTrackCuts::kBinLabel[] ={"TrkBit",
  45                                                    "TrkCuts",
  46                                                    "TrkEta",
  47                                                    "TrkDCA",
  48                                                    "TrkP",
  49                                                    "TrkPt",
  50                                                    "TrkPtTOF",
  51                                                    "TOFMatching",
  52                                                    "kTOFout",
  53                                                    "kTIME",
  54                                                    "kTOFpid",
  55                                                    "Accepted"};
  56
  57
  58 ClassImp(AliSpectraBothTrackCuts)
  59
  60
  61 AliSpectraBothTrackCuts::AliSpectraBothTrackCuts(const char *name) : TNamed(name, "AOD Track Cuts"), fIsSelected(0), fTrackBits(0), fMinTPCcls(0), fEtaCutMin(0), fEtaCutMax(0), fDCACut(0), fPCut(0), fPtCut(0), fYCutMax(0),fYCutMin(0),
  62   fPtCutTOFMatching(0),fAODtrack(kotherobject), fHashitinSPD1(0),fusedadditionalcuts(kTRUE),
  63 fHistoCuts(0), fHistoNSelectedPos(0), fHistoNSelectedNeg(0), fHistoNMatchedPos(0), fHistoNMatchedNeg(0), fHistoEtaPhiHighPt(0), fHistoNclustersITS(0),fTrack(0),fCuts(0)
  64
  65 {
  66   Bool_t oldStatus = TH1::AddDirectoryStatus();
  67   TH1::AddDirectory(kFALSE);
  68   // Constructor
  69   fHistoCuts = new TH1I("fTrkCuts", "Track Cuts", kNTrkCuts, -0.5, kNTrkCuts - 0.5);
  70   for(Int_t ibin=1;ibin<=kNTrkCuts;ibin++)fHistoCuts->GetXaxis()->SetBinLabel(ibin,kBinLabel[ibin-1]);
  71   //standard histo
  72   const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
  73   Int_t nbinsTempl=52;
  74
  75   fHistoNSelectedPos=new TH1F("fHistoNSelectedPos","fHistoNSelectedPos",nbinsTempl,templBins);
  76   fHistoNSelectedPos->GetXaxis()->SetTitle("P_{T} (GeV / c)");
  77   fHistoNSelectedNeg=new TH1F("fHistoNSelectedNeg","fHistoNSelectedNeg",nbinsTempl,templBins);
  78   fHistoNSelectedNeg->GetXaxis()->SetTitle("P_{T} (GeV / c)");
  79   fHistoNMatchedPos=new TH1F("fHistoNMatchedPos","fHistoNMatchedPos",nbinsTempl,templBins);
  80   fHistoNMatchedPos->GetXaxis()->SetTitle("P_{T} (GeV / c)");
  81   fHistoNMatchedNeg=new TH1F("fHistoNMatchedNeg","fHistoNMatchedNeg",nbinsTempl,templBins);
  82   fHistoNMatchedNeg->GetXaxis()->SetTitle("P_{T} (GeV / c)");
  83   fHistoEtaPhiHighPt=new TH2F("fHistoEtaPhiHighPt","fHistoEtaPhiHighPt",200,-1,1,400,0,7);
  84   fHistoEtaPhiHighPt->SetXTitle("eta");
  85   fHistoEtaPhiHighPt->SetYTitle("phi");
  86   fHistoNclustersITS=new TH1F("fHistoNclustersITS","fHistoNclustersITS;N;ITSLayer",6,-0.5,5.5);
  87
  88   fEtaCutMin = -100000.0; // default value of eta cut ~ no cut
  89   fEtaCutMax = 100000.0; // default value of eta cut ~ no cut
  90   fDCACut = 100000.0; // default value of dca cut ~ no cut
  91   fPCut = 100000.0; // default value of p cut ~ no cut
  92   fPtCut = 100000.0; // default value of pt cut ~ no cut
  93   fPtCutTOFMatching=0.6; //default value fot matching with TOF
  94   fYCutMax       = 100000.0; // default value of y cut ~ no cut
  95   fYCutMin       = -100000.0; // default value of y cut ~ no cut
  96   fMinTPCcls=70; // ncls in TPC
  97    TH1::AddDirectory(oldStatus);
  98
  99 }
 100 //_______________________________________________________
 101 AliSpectraBothTrackCuts::~AliSpectraBothTrackCuts()
 102 {
 103         if(fHistoCuts)
 104                 delete fHistoCuts;
 105         if(fHistoNSelectedPos)
 106                 delete fHistoNSelectedPos;
 107         if(fHistoNSelectedNeg)
 108                 delete fHistoNSelectedNeg;
 109         if(fHistoNMatchedPos)
 110                 delete fHistoNMatchedPos;
 111         if(fHistoNMatchedNeg)
 112                 delete fHistoNMatchedNeg;
 113         if(fHistoEtaPhiHighPt)
 114                 delete fHistoEtaPhiHighPt;
 115         if(fHistoNclustersITS)
 116                 delete fHistoNclustersITS;
 117
 118
 119
 120 }
 121 //_______________________________________________________
 122 Bool_t AliSpectraBothTrackCuts::IsSelected(AliVTrack * track,Bool_t FillHistStat)
 123 {
 124 // Returns true if Track Cuts are selected and applied
 125   if (!track)
 126     {
 127       printf("ERROR: Could not receive track");
 128       return kFALSE;
 129     }
 130     fTrack = track;
 131    TString nameoftrack(track->ClassName());
 132     if(!nameoftrack.CompareTo("AliESDtrack"))
 133                 fAODtrack=kESDobject;
 134         else if(!nameoftrack.CompareTo("AliAODTrack"))
 135                 fAODtrack=kAODobject;
 136         else
 137                 fAODtrack=kotherobject;
 138   if(!CheckTrackType()){
 139     return kFALSE;
 140   }
 141   if(FillHistStat)fHistoCuts->Fill(kTrkBit);
 142   if(!CheckTrackCuts()){
 143     return kFALSE;
 144   }
 145   if(FillHistStat)fHistoCuts->Fill(kTrkCuts);
 146   if(!CheckEtaCut()){
 147     return kFALSE;
 148   }
 149   if(FillHistStat)fHistoCuts->Fill(kTrkEta);
 150   if(!CheckDCACut()){
 151     return kFALSE;
 152   }
 153   if(FillHistStat)fHistoCuts->Fill(kTrkDCA);
 154   if(!CheckPCut()){
 155     return kFALSE;
 156   }
 157   if(FillHistStat)fHistoCuts->Fill(kTrkP);
 158   if(!CheckPtCut()){
 159     return kFALSE;
 160   }
 161   if(FillHistStat)fHistoCuts->Fill(kTrkPt);
 162   if(!CheckTOFMatching(FillHistStat)){
 163     return kFALSE;
 164   }
 165   if(FillHistStat)fHistoCuts->Fill(kAccepted);
 166   //Printf("-------- %d,%d",kTOFMatching,kAccepted);
 167
 168   return kTRUE;
 169 }
 170 //_________________________________________________________
 171
 172 Bool_t AliSpectraBothTrackCuts::CheckTrackType()
 173 {
 174   // Check track Type
 175   if(fAODtrack==kESDobject)
 176   {
 177         AliESDtrack* esdtrack=dynamic_cast<AliESDtrack*>(fTrack);
 178         if(!esdtrack)
 179                 return kFALSE;
 180         if(fCuts->AcceptTrack(esdtrack)) return kTRUE;
 181                 return kFALSE;
 182  }
 183   else if(fAODtrack==kAODobject)
 184   {
 185         AliAODTrack* aodtrack=dynamic_cast<AliAODTrack*>(fTrack);
 186         if(!aodtrack)
 187                 return kFALSE;
 188         if (aodtrack->TestFilterBit(fTrackBits)) return kTRUE;
 189                 return kFALSE;
 190   }
 191
 192   else
 193         return kFALSE;
 194
 195 }
 196 //_________________________________________________________
 197
 198 Bool_t AliSpectraBothTrackCuts::CheckTrackCuts()
 199 {
 200   // Check additional track Cuts
 201   Bool_t PassTrackCuts=kTRUE;
 202   if(!fusedadditionalcuts)
 203         return PassTrackCuts;
 204   AliAODTrack* aodtrack=0;
 205   AliESDtrack* esdtrack=0;
 206   if(fAODtrack==kESDobject)
 207   {
 208         esdtrack=dynamic_cast<AliESDtrack*>(fTrack);
 209         if(!esdtrack)
 210                 return kFALSE;
 211         if (!esdtrack->HasPointOnITSLayer(0) && !esdtrack->HasPointOnITSLayer(1))PassTrackCuts=kFALSE; //FIXME 1 SPD for the moment
 212         if (fHashitinSPD1&&!esdtrack->HasPointOnITSLayer(0)) PassTrackCuts=kFALSE;
 213         if (esdtrack->GetTPCNcls()<fMinTPCcls)PassTrackCuts=kFALSE;
 214         if(!esdtrack->IsOn(AliESDtrack::kTPCrefit))PassTrackCuts=kFALSE;
 215         if(!esdtrack->IsOn(AliESDtrack::kITSrefit))PassTrackCuts=kFALSE;
 216         if(PassTrackCuts)
 217         {
 218                 for(int i=0;i<6;i++)
 219                         if(esdtrack->HasPointOnITSLayer(i))
 220                                 fHistoNclustersITS->Fill(i);
 221         }
 222   }
 223   else if (fAODtrack==kAODobject)
 224         {
 225         aodtrack=dynamic_cast<AliAODTrack*>(fTrack);
 226         if(!aodtrack)
 227                 return kFALSE;
 228         if (!aodtrack->HasPointOnITSLayer(0) && !aodtrack->HasPointOnITSLayer(1))PassTrackCuts=kFALSE; //FIXME 1 SPD for the moment
 229         if (fHashitinSPD1&&!aodtrack->HasPointOnITSLayer(0)) PassTrackCuts=kFALSE;
 230         if (aodtrack->GetTPCNcls()<fMinTPCcls)PassTrackCuts=kFALSE;
 231         if(!aodtrack->IsOn(AliAODTrack::kTPCrefit))PassTrackCuts=kFALSE;
 232         if(!aodtrack->IsOn(AliAODTrack::kITSrefit))PassTrackCuts=kFALSE;
 233         if(PassTrackCuts)
 234         {
 235                 for(int i=0;i<6;i++)
 236                         if(aodtrack->HasPointOnITSLayer(i))
 237                                 fHistoNclustersITS->Fill(i);
 238         }
 239   }
 240   else
 241         return kFALSE;
 242
 243
 244
 245   return PassTrackCuts;
 246 }
 247 //________________________________________________________
 248 Bool_t AliSpectraBothTrackCuts::CheckEtaCut()
 249 {
 250    // Check eta cut
 251    if (fTrack->Eta() < fEtaCutMax && fTrack->Eta() > fEtaCutMin) return kTRUE;
 252     return kFALSE;
 253 }
 254
 255 Bool_t AliSpectraBothTrackCuts::CheckYCut(BothParticleSpecies_t species)
 256 {
 257   // check if the rapidity is within the set range
 258   Double_t y;
 259
 260   Double_t pz=fTrack->Pz();
 261   Double_t p=fTrack->P();
 262   Double_t mass=-1.0;
 263   /*
 264   if (species == kSpProton) { y = fTrack->Y(9.38271999999999995e-01); }
 265   if ( species == kSpKaon ) { y = fTrack->Y(4.93676999999999977e-01); }
 266   if ( species == kSpPion)  { y = fTrack->Y(1.39570000000000000e-01); }
 267
 268   */
 269     if (species == kSpProton) { mass=9.38271999999999995e-01; }
 270   if ( species == kSpKaon ) { mass=4.93676999999999977e-01; }
 271   if ( species == kSpPion)  { mass=1.39570000000000000e-01; }
 272   if(mass<0.0)
 273         y =-999.0 ;
 274   else
 275         y=0.5*TMath::Log((TMath::Sqrt(mass*mass+p*p)+pz)/(TMath::Sqrt(mass*mass+p*p)-pz));
 276   if (y > fYCutMax || y<fYCutMin||y < -998.) return kFALSE;
 277         return kTRUE;
 278 }
 279 //_______________________________________________________
 280 Bool_t AliSpectraBothTrackCuts::CheckDCACut()
 281 {
 282    // Check DCA cut
 283  // if (TMath::Abs(fTrack->DCA()) < fDCACut) return kTRUE; //FIXME for newest AOD fTrack->DCA() always gives -999
 284
 285     AliAODTrack* aodtrack=0;
 286   AliESDtrack* esdtrack=0;
 287   if(fAODtrack==kESDobject)
 288   {
 289         esdtrack=dynamic_cast<AliESDtrack*>(fTrack);
 290         if(!esdtrack)
 291                 return kFALSE;
 292         Float_t dcaxy=0.0;
 293         Float_t dcaz=0.0;
 294         esdtrack->GetImpactParameters(dcaxy,dcaz);
 295         if (TMath::Abs(dcaxy) < fDCACut)
 296                 return kTRUE;
 297         else
 298                 return kFALSE;
 299    }
 300   else if (fAODtrack==kAODobject)
 301   {
 302         aodtrack=dynamic_cast<AliAODTrack*>(fTrack);
 303         if(!aodtrack)
 304                 return kFALSE;
 305         if (TMath::Abs(aodtrack->DCA()) < fDCACut) return kTRUE;
 306         else
 307                 return kFALSE;
 308
 309    }
 310         else
 311         return kFALSE;
 312 }
 313 //________________________________________________________
 314 Bool_t AliSpectraBothTrackCuts::CheckPCut()
 315 {
 316    // Check P cut
 317    if (fTrack->P() < fPCut) return kTRUE;
 318    return kFALSE;
 319 }
 320 //_______________________________________________________
 321 Bool_t AliSpectraBothTrackCuts::CheckPtCut()
 322 {
 323     // check Pt cut
 324 //    if ((fTrack->Pt() < fPtCut) && (fTrack->Pt() > 0.3 )) return kTRUE;
 325    if (fTrack->Pt() < fPtCut) return kTRUE;
 326     return kFALSE;
 327 }
 328
 329 //_______________________________________________________
 330 Bool_t AliSpectraBothTrackCuts::CheckTOFMatching(Bool_t FillHistStat)
 331 {
 332   // check Pt cut
 333   //    if ((fTrack->Pt() < fPtCut) && (fTrack->Pt() > 0.3 )) return kTRUE;
 334
 335   if (fTrack->Pt() < fPtCutTOFMatching) return kTRUE;
 336   else{
 337     if(FillHistStat)fHistoCuts->Fill(kTrkPtTOF);
 338     if(fTrack->Charge()>0)fHistoNSelectedPos->Fill(fTrack->Pt());
 339     else fHistoNSelectedNeg->Fill(fTrack->Pt());
 340     UInt_t status;
 341     status=fTrack->GetStatus();
 342     if((status&AliAODTrack::kTOFout)&&FillHistStat)fHistoCuts->Fill(kTrTOFout);
 343     if((status&AliAODTrack::kTIME)&&FillHistStat)fHistoCuts->Fill(kTrTIME);
 344     if((status&AliAODTrack::kTOFpid)&&FillHistStat)fHistoCuts->Fill(kTrTOFpid);
 345
 346     if((status&AliAODTrack::kTOFout)==0 || (status&AliAODTrack::kTIME)==0){//kTOFout and kTIME
 347       return kFALSE;
 348     }
 349     if(FillHistStat)fHistoCuts->Fill(kTOFMatching);
 350     if(fTrack->Charge()>0)fHistoNMatchedPos->Fill(fTrack->Pt());
 351     else fHistoNMatchedNeg->Fill(fTrack->Pt());
 352     if(fTrack->Pt()>1.5){
 353       //fHistoEtaPhiHighPt->Fill(fTrack->GetOuterParam()->Eta(),fTrack->GetOuterParam()->Phi());
 354       //Printf("AliExternalTrackParam * extpar=(AliExternalTrackParam*)fTrack->GetOuterParam();");
 355       //AliExternalTrackParam * extpar=(AliExternalTrackParam*)fTrack->GetOuterParam();
 356       fHistoEtaPhiHighPt->Fill(fTrack->Eta(),fTrack->Phi());
 357       //Printf("fHistoEtaPhiHighPt->Fill(extpar->Eta(),extpar->Phi());");
 358       //fHistoEtaPhiHighPt->Fill(extpar->Eta(),extpar->Phi());
 359       //delete extpar;
 360     }
 361     return kTRUE;
 362   }
 363 }
 364 //_______________________________________________________
 365 void AliSpectraBothTrackCuts::PrintCuts() const
 366 {
 367   // Print cuts
 368     cout << "Track Cuts" << endl;
 369     cout << " > TrackBit\t" << fTrackBits << endl;
 370     cout << " > Eta cut\t" << fEtaCutMin <<","<< fEtaCutMax << endl;
 371     cout << " > DCA cut\t" << fDCACut << endl;
 372     cout << " > P cut\t" << fPCut << endl;
 373     cout << " > Pt cut \t" << fPtCut << endl;
 374     cout << " > TPC cls \t" << fMinTPCcls << endl;
 375 }
 376 //_______________________________________________________
 377 void AliSpectraBothTrackCuts::SetTrackType(UInt_t bit)
 378 {
 379    // Set the type of track to be used. The argument should be the bit number. The mask is produced automatically.
 380    fTrackBits = (0x1 << (bit - 1));
 381 }
 382 //_______________________________________________________
 383
 384 Long64_t AliSpectraBothTrackCuts::Merge(TCollection* list)
 385 {
 386   // Merge a list of AliSpectraBothTrackCuts objects with this.
 387   // Returns the number of merged objects (including this).
 388
 389   //  AliInfo("Merging");
 390
 391
 392   if (!list)
 393     return 0;
 394
 395   if (list->IsEmpty())
 396     return 1;
 397
 398   TIterator* iter = list->MakeIterator();
 399   TObject* obj;
 400
 401   // collections of all histograms
 402   TList collections;//FIXME we should only 1 collection
 403   TList collections_histoNSelectedPos;
 404   TList collections_histoNSelectedNeg;
 405   TList collections_histoNMatchedPos;
 406   TList collections_histoNMatchedNeg;
 407   TList collections_histoEtaPhiHighPt;
 408
 409   Int_t count = 0;
 410
 411   while ((obj = iter->Next())) {
 412     AliSpectraBothTrackCuts* entry = dynamic_cast<AliSpectraBothTrackCuts*> (obj);
 413     if (entry == 0)
 414       continue;
 415
 416     TH1I * histo = entry->GetHistoCuts();
 417     collections.Add(histo);
 418     TH1F * histoNSelectedPos = entry->GetHistoNSelectedPos();
 419     collections_histoNSelectedPos.Add(histoNSelectedPos);
 420     TH1F * histoNSelectedNeg = entry->GetHistoNSelectedNeg();
 421     collections_histoNSelectedNeg.Add(histoNSelectedNeg);
 422     TH1F * histoNMatchedPos = entry->GetHistoNMatchedPos();
 423     collections_histoNMatchedPos.Add(histoNMatchedPos);
 424     TH1F * histoNMatchedNeg = entry->GetHistoNMatchedNeg();
 425     collections_histoNMatchedNeg.Add(histoNMatchedNeg);
 426     TH2F * histoEtaPhiHighPt = entry->GetHistoEtaPhiHighPt();
 427     collections_histoEtaPhiHighPt.Add(histoEtaPhiHighPt);
 428     count++;
 429   }
 430
 431   fHistoCuts->Merge(&collections);
 432   fHistoNSelectedPos->Merge(&collections_histoNSelectedPos);
 433   fHistoNSelectedNeg->Merge(&collections_histoNSelectedNeg);
 434   fHistoNMatchedPos->Merge(&collections_histoNMatchedPos);
 435   fHistoNMatchedNeg->Merge(&collections_histoNMatchedNeg);
 436   fHistoEtaPhiHighPt->Merge(&collections_histoEtaPhiHighPt);
 437
 438   delete iter;
 439
 440   return count+1;
 441 }
 442