]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG2/SPECTRA/AliProtonSpectraCorrection.cxx
Request by Martin: added flag for big output
[u/mrichter/AliRoot.git] / 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 }