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",22);
 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",17);
 119   strncpy(nameCF,"MC particle generation cuts",28);
 120   AliCFParticleGenCuts* mcGenCutsProtons = new AliCFParticleGenCuts(titleCF,nameCF);
 121   mcGenCutsProtons->SetRequireIsPrimary();
 122   mcGenCutsProtons->SetRequirePdgCode(2212);
 123
 124   strncpy(titleCF,"mcGenCutsAntiProtons",21);
 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",10);
 138   strncpy(nameCF,"Acceptance cuts",16);
 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",19);
 148   strncpy(nameCF,"rec-level kine cuts",20);
 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",23);
 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 (mcPart){
 210       if (fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 211                                                mcPart)) {
 212         containerInput[0] = (Float_t)mcPart->Eta();
 213         containerInput[1] = (Float_t)mcPart->Pt();
 214         //fill the container for Gen-level selection
 215         fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 216                                                         kStepGenerated);
 217         //check the Acceptance-level cuts
 218         if (fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,
 219                                                  mcPart)) {
 220           //fill the container for Acceptance-level selection
 221           fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 222                                                           kStepReconstructible);
 223         }//acceptance cuts - protons
 224       }//MC level cuts - protons
 225
 226       //Antiprotons
 227       //check the MC-level cuts
 228       if (fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 229                                                    mcPart)) {
 230         containerInput[0] = (Float_t)mcPart->Eta();
 231         containerInput[1] = (Float_t)mcPart->Pt();
 232         //fill the container for Gen-level selection
 233         fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 234                                                             kStepGenerated);
 235         //check the Acceptance-level cuts
 236         if (fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,
 237                                                      mcPart)) {
 238           //fill the container for Acceptance-level selection
 239           fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 240                                                               kStepReconstructible);
 241         }//acceptance cuts - antiprotons
 242       }//MC level cuts - antiprotons
 243     }//MC particle condition
 244   }//loop over MC particles
 245
 246
 247   //__________________________________________________________//
 248   //ESD track loop
 249   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 250     AliESDtrack *track = dynamic_cast<AliESDtrack *>(esd->GetTrack(iTrack));
 251     if(!track) continue;
 252
 253     // is track associated to particle ?
 254     Int_t label = track->GetLabel();
 255     if (label<0) continue;
 256     AliMCParticle *mcPart  = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
 257     if (!mcPart) continue;
 258
 259     if((fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC)||(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kHybrid)) {
 260       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 261       if(!tpcTrack) continue;
 262     }//Hybrid or standalone TPC
 263
 264     //Protons
 265     // check if this track was part of the signal - primaries
 266     if (fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 267                                              mcPart)) {
 268       //fill the container - reconstructed protons
 269       containerInput[0] = mcPart->Eta();
 270       containerInput[1] = mcPart->Pt();
 271       fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 272                                                       kStepReconstructed);
 273       //fill the container - identified protons
 274       if(fProtonAnalysisBase->IsProton(track)) {
 275         containerInput[0] = mcPart->Eta();
 276         containerInput[1] = mcPart->Pt();
 277         fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 278                                                         kStepIdentified);
 279
 280         //fill the container - survived protons
 281         //track cuts
 282         if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track)) continue;
 283         if(!fProtonAnalysisBase->IsAccepted(track)) continue;
 284         //track outside the analyzed y-Pt
 285         if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue;
 286         containerInput[0] = mcPart->Eta();
 287         containerInput[1] = mcPart->Pt();
 288         fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 289                                                         kStepSelected);
 290       }
 291     }//MC primaries - protons
 292
 293     //Antiprotons
 294     // check if this track was part of the signal - primaries
 295     if (fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 296                                                  mcPart)) {
 297       //fill the container - reconstructed antiprotons
 298       containerInput[0] = mcPart->Eta();
 299       containerInput[1] = mcPart->Pt();
 300       fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 301                                                           kStepReconstructed);
 302
 303       //fill the container - identified antiprotons
 304       if(fProtonAnalysisBase->IsProton(track)) {
 305         containerInput[0] = mcPart->Eta();
 306         containerInput[1] = mcPart->Pt();
 307         fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 308                                                             kStepIdentified);
 309
 310         //fill the container - survived antiprotons
 311         //track cuts
 312         if(!fProtonAnalysisBase->IsPrimary(esd,vertex,track)) continue;
 313         if(!fProtonAnalysisBase->IsAccepted(track)) continue;
 314         //track outside the analyzed y-Pt
 315         if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue;
 316         containerInput[0] = mcPart->Eta();
 317         containerInput[1] = mcPart->Pt();
 318         fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 319                                                             kStepSelected);
 320       }
 321     }//MC primaries - antiprotons
 322   }//track loop
 323
 324   //if(fProtonAnalysisBase->GetDebugMode())
 325   //Printf("Initial number of tracks: %d | Identified (anti)protons: %d - %d | Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);
 326 }
 327
 328 //_________________________________________________________________________//
 329 void AliProtonSpectraCorrection::FillCorrectionMaps(AliAODEvent *fAOD) {
 330   // Analysis from AOD: to be implemented...
 331   // in the near future.
 332   fAOD->Print();
 333 }