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