PWG2/SPECTRA/AliProtonAbsorptionCorrection.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 "AliProtonAbsorptionCorrection.h"
  30 #include "AliProtonAnalysisBase.h"
  31
  32 ClassImp(AliProtonAbsorptionCorrection)
  33
  34 //____________________________________________________________________//
  35 AliProtonAbsorptionCorrection::AliProtonAbsorptionCorrection() :
  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 AliProtonAbsorptionCorrection::~AliProtonAbsorptionCorrection() {
  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 AliProtonAbsorptionCorrection::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",27);
 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 AliProtonAbsorptionCorrection::FillAbsorptionMaps(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     if(!mcPart) continue;
 207
 208     Double_t vz = mcPart->Zv();
 209     if (TMath::Abs(vz) > 50.) continue;//exclude particles generated out of the acceptance
 210     if(TMath::Abs(mcPart->Eta()) > 1.0) continue;
 211
 212     //Protons
 213     //check the MC-level cuts
 214     if (fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 215                                              mcPart)) {
 216       if(fProtonAnalysisBase->GetEtaMode())
 217         containerInput[0] = (Float_t)mcPart->Eta();
 218       else
 219         containerInput[0] = (Float_t)fProtonAnalysisBase->Rapidity(mcPart->Px(),mcPart->Pz(),mcPart->Pz());
 220       containerInput[1] = (Float_t)mcPart->Pt();
 221       //fill the container for Gen-level selection
 222       fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 223                                                       kStepGenerated);
 224       //check the Acceptance-level cuts
 225       if (fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,
 226                                                mcPart)) {
 227         //fill the container for Acceptance-level selection
 228         fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 229                                                         kStepReconstructible);
 230       }//acceptance cuts - protons
 231     }//MC level cuts - protons
 232
 233     //Antiprotons
 234     //check the MC-level cuts
 235     if (fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 236                                                  mcPart)) {
 237       if(fProtonAnalysisBase->GetEtaMode())
 238         containerInput[0] = (Float_t)mcPart->Eta();
 239       else
 240         containerInput[0] = (Float_t)fProtonAnalysisBase->Rapidity(mcPart->Px(),mcPart->Pz(),mcPart->Pz());
 241       containerInput[1] = (Float_t)mcPart->Pt();
 242       //fill the container for Gen-level selection
 243       fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 244                                                           kStepGenerated);
 245       //check the Acceptance-level cuts
 246       if (fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,
 247                                                    mcPart)) {
 248         //fill the container for Acceptance-level selection
 249         fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 250                                                             kStepReconstructible);
 251       }//acceptance cuts - antiprotons
 252     }//MC level cuts - antiprotons
 253   }//loop over MC particles
 254
 255
 256   //__________________________________________________________//
 257   //ESD track loop
 258   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 259     AliESDtrack *track = dynamic_cast<AliESDtrack *>(esd->GetTrack(iTrack));
 260     if(!track) continue;
 261
 262     // is track associated to particle ?
 263     Int_t label = track->GetLabel();
 264     if (label < 0) continue;
 265     if(label > mcEvent->GetNumberOfTracks()) continue;
 266
 267     AliMCParticle *mcPart  = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
 268     if(!mcPart) continue;
 269
 270     Double_t vz = mcPart->Zv();
 271     if (TMath::Abs(vz) > 50.) continue;//exclude particles generated out of the acceptance
 272     if(TMath::Abs(mcPart->Eta()) > 1.0) continue;
 273
 274     if((fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC)||(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kHybrid)) {
 275       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 276       if(!tpcTrack) continue;
 277     }//Hybrid or standalone TPC
 278
 279     //Protons
 280     // check if this track was part of the signal
 281     if (fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 282                                              mcPart)) {
 283       //fill the container
 284       if(fProtonAnalysisBase->GetEtaMode())
 285         containerInput[0] = mcPart->Eta();
 286       else
 287         containerInput[0] = (Float_t)fProtonAnalysisBase->Rapidity(mcPart->Px(),mcPart->Pz(),mcPart->Pz());
 288       containerInput[1] = mcPart->Pt();
 289       fCFManagerProtons->GetParticleContainer()->Fill(containerInput,
 290                                                       kStepReconstructed) ;
 291     }//MC primaries - protons
 292
 293     //Antiprotons
 294     // check if this track was part of the signal
 295     if (fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,
 296                                                  mcPart)) {
 297       //fill the container
 298       if(fProtonAnalysisBase->GetEtaMode())
 299         containerInput[0] = mcPart->Eta();
 300       else
 301         containerInput[0] = (Float_t)fProtonAnalysisBase->Rapidity(mcPart->Px(),mcPart->Pz(),mcPart->Pz());
 302       containerInput[1] = mcPart->Pt();
 303       fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,
 304                                                           kStepReconstructed);
 305     }//MC primaries - antiprotons
 306   }//track loop
 307
 308   //if(fProtonAnalysisBase->GetDebugMode())
 309   //Printf("Initial number of tracks: %d | Identified (anti)protons: %d - %d | Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);
 310 }
 311
 312 //_________________________________________________________________________//
 313 void AliProtonAbsorptionCorrection::FillAbsorptionMaps(AliAODEvent *fAOD) {
 314   // Analysis from AOD: to be implemented...
 315   // in the near future.
 316   fAOD->Print();
 317 }