PWG2/SPECTRA/AliProtonSpectraCorrection.cxx

   1 #include <Riostream.h>
   2 #include <TFile.h>
   3 #include <TSystem.h>
   4 #include <TF1.h>
   5 #include <TH2D.h>
   6 #include <TH1D.h>
   7 #include <TH1I.h>
   8 #include <TParticle.h>
   9
  10 #include <AliExternalTrackParam.h>
  11 #include <AliAODEvent.h>
  12 #include <AliESDEvent.h>
  13 //#include <AliLog.h>
  14 #include <AliPID.h>
  15 #include <AliStack.h>
  16 #include <AliCFContainer.h>
  17 #include <AliCFEffGrid.h>
  18 #include <AliCFDataGrid.h>
  19 #include <AliCFManager.h>
  20 #include <AliCFTrackKineCuts.h>
  21 #include <AliCFParticleGenCuts.h>
  22 #include <AliCFAcceptanceCuts.h>
  23 #include <AliMCEvent.h>
  24 //#include <AliESDVertex.h>
  25 class AliLog;
  26 class AliESDVertex;
  27
  28 #include "AliProtonSpectraCorrection.h"
  29 #include "AliProtonAnalysisBase.h"
  30
  31 ClassImp(AliProtonSpectraCorrection)
  32
  33 //____________________________________________________________________//
  34 AliProtonSpectraCorrection::AliProtonSpectraCorrection() :
  35   TObject(), fProtonAnalysisBase(0),
  36   fNBinsY(0), fMinY(0), fMaxY(0),
  37   fNBinsPt(0), fMinPt(0), fMaxPt(0),
  38   fCFManagerProtons(0), fCFManagerAntiProtons(0) {
  39   //fProtonContainer(0), fAntiProtonContainer(0) {
  40   //Default constructor
  41 }
  42
  43 //____________________________________________________________________//
  44 AliProtonSpectraCorrection::~AliProtonSpectraCorrection() {
  45   //Default destructor
  46   if(fCFManagerProtons) delete fCFManagerProtons;
  47   if(fCFManagerAntiProtons) delete fCFManagerAntiProtons;
  48   if(fProtonAnalysisBase) delete fProtonAnalysisBase;
  49   //if(fProtonContainer) delete fProtonContainer;
  50   //if(fAntiProtonContainer) delete fAntiProtonContainer;
  51 }
  52
  53 //____________________________________________________________________//
  54 void AliProtonSpectraCorrection::InitAnalysisHistograms(Int_t nbinsY,
  55                                                            Float_t fLowY,
  56                                                            Float_t fHighY,
  57                                                            Int_t nbinsPt,
  58                                                            Float_t fLowPt,
  59                                                            Float_t fHighPt) {
  60   //Initializes the histograms
  61   fNBinsY = nbinsY;
  62   fMinY = fLowY;
  63   fMaxY = fHighY;
  64   fNBinsPt = nbinsPt;
  65   fMinPt = fLowPt;
  66   fMaxPt = fHighPt;
  67   const Int_t    mintrackrefsTPC = 1;
  68
  69   //=========================================================//
  70   //setting up the containers
  71   Int_t iBin[2];
  72   iBin[0] = nbinsY;
  73   iBin[1] = nbinsPt;
  74   Double_t *binLimY = new Double_t[nbinsY+1];
  75   Double_t *binLimPt = new Double_t[nbinsPt+1];
  76   //values for bin lower bounds
  77   for(Int_t i = 0; i <= nbinsY; i++)
  78     binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
  79   for(Int_t i = 0; i <= nbinsPt; i++)
  80     binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
  81
  82   //Proton container
  83   AliCFContainer *containerProtons = new AliCFContainer("containerProtons",
  84                                                         "container for protons",
  85                                                         kNSteps,2,iBin);
  86   containerProtons->SetBinLimits(0,binLimY); //rapidity or eta
  87   containerProtons->SetBinLimits(1,binLimPt); //pT
  88
  89   //Anti-proton container
  90   AliCFContainer *containerAntiProtons = new AliCFContainer("containerAntiProtons",
  91                                                             "container for antiprotons",
  92                                                             kNSteps,2,iBin);
  93   containerAntiProtons->SetBinLimits(0,binLimY); //rapidity or eta
  94   containerAntiProtons->SetBinLimits(1,binLimPt); //pT
  95
  96   //=========================================================//
  97   //Setting up the criteria for the generated particles
  98   AliCFTrackKineCuts *mcKineCutsProtons = new AliCFTrackKineCuts("mcKineCutsProtons",
  99                                                                  "MC-level kinematic cuts");
 100   mcKineCutsProtons->SetPtRange(fMinPt,fMaxPt);
 101   if(fProtonAnalysisBase->GetEtaMode())
 102     mcKineCutsProtons->SetEtaRange(fMinY,fMaxY);
 103   else
 104     mcKineCutsProtons->SetRapidityRange(fMinY,fMaxY);
 105   mcKineCutsProtons->SetChargeMC(1.0);
 106
 107   AliCFTrackKineCuts *mcKineCutsAntiProtons = new AliCFTrackKineCuts("mcKineCutsAntiProtons",
 108                                                                      "MC-level kinematic cuts");
 109   mcKineCutsAntiProtons->SetPtRange(fMinPt,fMaxPt);
 110   if(fProtonAnalysisBase->GetEtaMode())
 111     mcKineCutsAntiProtons->SetEtaRange(fMinY,fMaxY);
 112   else
 113     mcKineCutsAntiProtons->SetRapidityRange(fMinY,fMaxY);
 114   mcKineCutsAntiProtons->SetChargeMC(-1.0);
 115
 116   AliCFParticleGenCuts* mcGenCutsProtons = new AliCFParticleGenCuts("mcGenCutsProtons",
 117                                                                     "MC particle generation cuts");
 118   mcGenCutsProtons->SetRequireIsPrimary();
 119   mcGenCutsProtons->SetRequirePdgCode(2212);
 120   AliCFParticleGenCuts* mcGenCutsAntiProtons = new AliCFParticleGenCuts("mcGenCutsAntiProtons",
 121                                                                         "MC particle generation cuts");
 122   mcGenCutsAntiProtons->SetRequireIsPrimary();
 123   mcGenCutsAntiProtons->SetRequirePdgCode(-2212);
 124
 125   TObjArray* mcListProtons = new TObjArray(0);
 126   mcListProtons->AddLast(mcKineCutsProtons);
 127   mcListProtons->AddLast(mcGenCutsProtons);
 128   TObjArray* mcListAntiProtons = new TObjArray(0);
 129   mcListAntiProtons->AddLast(mcKineCutsAntiProtons);
 130   mcListAntiProtons->AddLast(mcGenCutsAntiProtons);
 131
 132   //Setting up the criteria for the reconstructible particles
 133   AliCFAcceptanceCuts *mcAccCuts = new AliCFAcceptanceCuts("mcAccCuts",
 134                                                            "Acceptance cuts");
 135   mcAccCuts->SetMinNHitTPC(mintrackrefsTPC);
 136   TObjArray* accList = new TObjArray(0);
 137   accList->AddLast(mcAccCuts);
 138
 139   //____________________________________________//
 140   //Setting up the criteria for the reconstructed tracks
 141   //____________________________________________//
 142   AliCFTrackKineCuts *recKineCutsProtons = new AliCFTrackKineCuts("recKineCutsProtons",
 143                                                                   "rec-level kine cuts");
 144   recKineCutsProtons->SetPtRange(fMinPt,fMaxPt);
 145   if(fProtonAnalysisBase->GetEtaMode())
 146     recKineCutsProtons->SetEtaRange(fMinY,fMaxY);
 147   else
 148     recKineCutsProtons->SetRapidityRange(fMinY,fMaxY);
 149   recKineCutsProtons->SetChargeRec(1.0);
 150
 151   //____________________________________________//
 152   AliCFTrackKineCuts *recKineCutsAntiProtons = new AliCFTrackKineCuts("recKineCutsAntiProtons",
 153                                                                       "rec-level kine cuts");
 154   recKineCutsAntiProtons->SetPtRange(fMinPt,fMaxPt);
 155   if(fProtonAnalysisBase->GetEtaMode())
 156     recKineCutsAntiProtons->SetEtaRange(fMinY,fMinY);
 157   else
 158     recKineCutsAntiProtons->SetRapidityRange(fMinY,fMaxY);
 159   recKineCutsAntiProtons->SetChargeRec(-1.0);
 160
 161   //____________________________________________//
 162   TObjArray* recListProtons = new TObjArray(0);
 163   recListProtons->AddLast(recKineCutsProtons);
 164   recListProtons->AddLast(mcGenCutsProtons);
 165
 166   TObjArray* recListAntiProtons = new TObjArray(0);
 167   recListAntiProtons->AddLast(recKineCutsAntiProtons);
 168   recListAntiProtons->AddLast(mcGenCutsAntiProtons);
 169
 170   //=========================================================//
 171   //CF manager - Protons
 172   fCFManagerProtons = new AliCFManager();
 173   fCFManagerProtons->SetParticleContainer(containerProtons);
 174   fCFManagerProtons->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListProtons);
 175   fCFManagerProtons->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
 176   fCFManagerProtons->SetParticleCutsList(AliCFManager::kPartRecCuts,recListProtons);
 177
 178   //CF manager - Protons
 179   fCFManagerAntiProtons = new AliCFManager();
 180   fCFManagerAntiProtons->SetParticleContainer(containerAntiProtons);
 181   fCFManagerAntiProtons->SetParticleCutsList(AliCFManager::kPartGenCuts,mcListAntiProtons);
 182   fCFManagerAntiProtons->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
 183   fCFManagerAntiProtons->SetParticleCutsList(AliCFManager::kPartRecCuts,recListAntiProtons);
 184 }
 185
 186 //_________________________________________________________________________//
 187 void AliProtonSpectraCorrection::FillCorrectionMaps(AliESDEvent *esd,
 188                                                     const AliESDVertex *vertex,
 189                                                     AliMCEvent *mcEvent) {
 190   //Function to fill the correction containers
 191   fCFManagerProtons->SetEventInfo(mcEvent);
 192   fCFManagerAntiProtons->SetEventInfo(mcEvent);
 193
 194   //Dummy container
 195   Double_t containerInput[2];
 196   //__________________________________________________________//
 197   //loop on the MC event
 198   for (Int_t ipart = 0; ipart < mcEvent->GetNumberOfTracks(); ipart++) {
 199     AliMCParticle *mcPart  = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(ipart));
 200
 201     //Protons
 202     //check the MC-level cuts
 203     if (fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 204                                              mcPart)) {
 205       containerInput[0] = (Float_t)mcPart->Eta();
 206       containerInput[1] = (Float_t)mcPart->Pt();
 207       //fill the container for Gen-level selection
 208       fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 209                                                       kStepGenerated);
 210       //check the Acceptance-level cuts
 211       if (fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,
 212                                                mcPart)) {
 213         //fill the container for Acceptance-level selection
 214         fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 215                                                         kStepReconstructible);
 216       }//acceptance cuts - protons
 217     }//MC level cuts - protons
 218
 219     //Antiprotons
 220     //check the MC-level cuts
 221     if (fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 222                                                  mcPart)) {
 223       containerInput[0] = (Float_t)mcPart->Eta();
 224       containerInput[1] = (Float_t)mcPart->Pt();
 225       //fill the container for Gen-level selection
 226       fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 227                                                           kStepGenerated);
 228       //check the Acceptance-level cuts
 229       if (fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,
 230                                                    mcPart)) {
 231         //fill the container for Acceptance-level selection
 232         fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 233                                                             kStepReconstructible);
 234       }//acceptance cuts - antiprotons
 235     }//MC level cuts - antiprotons
 236   }//loop over MC particles
 237
 238
 239   //__________________________________________________________//
 240   //ESD track loop
 241   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 242     AliESDtrack *track = dynamic_cast<AliESDtrack *>(esd->GetTrack(iTrack));
 243     if(!track) continue;
 244
 245     // is track associated to particle ?
 246     Int_t label = track->GetLabel();
 247     if (label<0) continue;
 248     AliMCParticle *mcPart  = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
 249
 250     if((fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC)||(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kHybrid)) {
 251       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 252       if(!tpcTrack) continue;
 253     }//Hybrid or standalone TPC
 254
 255     //Protons
 256     // check if this track was part of the signal - primaries
 257     if (fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 258                                              mcPart)) {
 259       //fill the container - reconstructed protons
 260       containerInput[0] = mcPart->Eta();
 261       containerInput[1] = mcPart->Pt();
 262       fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 263                                                       kStepReconstructed);
 264       //fill the container - identified protons
 265       if(fProtonAnalysisBase->IsProton(track)) {
 266         containerInput[0] = mcPart->Eta();
 267         containerInput[1] = mcPart->Pt();
 268         fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 269                                                         kStepIdentified);
 270
 271         //fill the container - survived protons
 272         //track cuts
 273         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) continue;
 274         //track outside the analyzed y-Pt
 275         if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue;
 276         containerInput[0] = mcPart->Eta();
 277         containerInput[1] = mcPart->Pt();
 278         fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 279                                                         kStepSelected);
 280       }
 281     }//MC primaries - protons
 282
 283     //Antiprotons
 284     // check if this track was part of the signal - primaries
 285     if (fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 286                                                  mcPart)) {
 287       //fill the container - reconstructed antiprotons
 288       containerInput[0] = mcPart->Eta();
 289       containerInput[1] = mcPart->Pt();
 290       fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 291                                                           kStepReconstructed);
 292
 293       //fill the container - identified antiprotons
 294       if(fProtonAnalysisBase->IsProton(track)) {
 295         containerInput[0] = mcPart->Eta();
 296         containerInput[1] = mcPart->Pt();
 297         fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 298                                                             kStepIdentified);
 299
 300         //fill the container - survived antiprotons
 301         //track cuts
 302         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) continue;
 303         //track outside the analyzed y-Pt
 304         if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue;
 305         containerInput[0] = mcPart->Eta();
 306         containerInput[1] = mcPart->Pt();
 307         fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 308                                                             kStepSelected);
 309       }
 310     }//MC primaries - antiprotons
 311   }//track loop
 312
 313   //if(fProtonAnalysisBase->GetDebugMode())
 314   //Printf("Initial number of tracks: %d | Identified (anti)protons: %d - %d | Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);
 315 }
 316
 317 //_________________________________________________________________________//
 318 void AliProtonSpectraCorrection::FillCorrectionMaps(AliAODEvent *fAOD) {
 319   // Analysis from AOD: to be implemented...
 320   // in the near future.
 321   fAOD->Print();
 322 }