Update PR task: drathee
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidRatioDCA.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: ALICE Offline.                                                 *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 //=========================================================================//
18 //             AliEbyE Analysis for Particle Ratio Fluctuation             //
19 //                   Deepika Rathee  | Satyajit Jena                       //
20 //                   drathee@cern.ch | sjena@cern.ch                       //
21 //                  Date: Wed Jul  9 18:38:30 CEST 2014                    //
22 //          New approch to find particle ratio to reduce memory            //
23 //                             (Test Only)                                 //
24 //=========================================================================//
25
26 #include "TMath.h"
27 #include "TAxis.h"
28
29 #include "AliESDEvent.h"
30 #include "AliStack.h"
31 #include "AliMCEvent.h"
32 #include "AliESDtrackCuts.h"
33
34 #include "AliAODEvent.h"
35 #include "AliAODMCParticle.h"
36
37 #include "AliEbyEPidRatioDCA.h"
38
39 using namespace std;
40
41 ClassImp(AliEbyEPidRatioDCA)
42
43 //________________________________________________________________________
44 AliEbyEPidRatioDCA::AliEbyEPidRatioDCA() :
45   AliEbyEPidRatioBase("DCA", "DCA"),
46   fESDTrackCutsBkg(NULL),
47   fHnDCA(NULL) {
48   // Constructor   
49
50   AliLog::SetClassDebugLevel("AliEbyEPidRatioDCA",10);
51 }
52
53 //________________________________________________________________________
54 AliEbyEPidRatioDCA::~AliEbyEPidRatioDCA() {
55   // Destructor
56 }
57
58 //________________________________________________________________________
59 void AliEbyEPidRatioDCA::Process() {
60    Float_t etaRange[2];
61   fESDTrackCutsBkg->GetEtaRange(etaRange[0],etaRange[1]);
62
63   Float_t ptRange[2];
64   fESDTrackCutsBkg->GetPtRange(ptRange[0],ptRange[1]); 
65
66   // -- Track Loop
67   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
68     
69     AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack)); 
70
71     // -- Check if track is accepted for basic parameters
72     if (!fHelper->IsTrackAcceptedBasicCharged(track))
73       continue;
74     
75     // -- Check if accepted - ESD
76     if (fESD && !fESDTrackCutsBkg->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
77       continue;
78     
79     // -- Check if accepted - AOD
80     if (fAOD){
81       AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
82       
83       if (!trackAOD) {
84         AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
85         continue;
86       }
87       if (!trackAOD->TestFilterBit(fAODtrackCutBit))
88         continue;
89
90       // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
91       if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
92         continue;
93     }
94
95     Int_t gPdgCode = 0;
96     Int_t iPid = 0;
97     Double_t pid[3];
98     if      (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion)))  {  iPid = 1; gPdgCode = 211;}
99     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon)))  {  iPid = 2; gPdgCode = 321;}
100     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))){  iPid = 3; gPdgCode = 2212;}
101     else iPid = 0;
102
103     //  cout << " --- DCA ---- " << iPid << "  " << gPdgCode << endl;
104     Double_t yP;
105     if (!fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
106       continue;
107   
108     Bool_t isDCArAccepted = fHelper->IsTrackAcceptedDCA(track);
109
110     // -- Check if accepted with thighter DCA cuts
111     // ?!?!? How to mimic this in AODs?
112     if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
113       isDCArAccepted = kFALSE;
114     
115     // -- Check for contamination 
116     Int_t contIdx = (fIsMC) ? GetContIdxTrack(TMath::Abs(track->GetLabel()), track->Charge(), gPdgCode) : 1;
117     
118     // -- Get DCAs (dca_r, dca_z, sigma_xy, sigma_xy_z, sigma_z)
119     Float_t dca[2], cov[3]; // 
120     if (fESD)
121       (dynamic_cast<AliESDtrack*>(track))->GetImpactParameters(dca, cov);
122     else 
123       dca[0] = 1.; 
124
125     // -- Fill THnSparse 
126     
127     if(iPid != 0) {   
128       Double_t hnDCA[10] = {fCentralityBin,0, 
129                             track->Eta(), 
130                             yP, 
131                             track->Phi(), 
132                             track->Pt(), 
133                             track->Charge(), 
134                             contIdx,
135                             isDCArAccepted, 
136                             dca[0]};
137       fHnDCA->Fill(hnDCA);
138     }      
139     
140     Double_t hnDCA[10] = {fCentralityBin, iPid,
141                           track->Eta(), 
142                           yP, 
143                           track->Phi(), 
144                           track->Pt(), 
145                           track->Charge(), 
146                           contIdx,
147                           isDCArAccepted, 
148                           dca[0]};
149       fHnDCA->Fill(hnDCA);
150
151
152     } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
153
154   return;
155 }      
156
157 //________________________________________________________________________
158 void AliEbyEPidRatioDCA::CreateHistograms() {
159   Int_t    binHnDCA[11] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent,4,
160                            AliEbyEPidRatioHelper::fgkfHistNBinsEta,       
161                            AliEbyEPidRatioHelper::fgkfHistNBinsRap,  
162                            AliEbyEPidRatioHelper::fgkfHistNBinsPhi,        
163                            AliEbyEPidRatioHelper::fgkfHistNBinsPt,   
164                            AliEbyEPidRatioHelper::fgkfHistNBinsSign,      
165                            4,  2,  77};      
166   
167   Double_t minHnDCA[11] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0],-0.5, 
168                            AliEbyEPidRatioHelper::fgkfHistRangeEta[0], 
169                            AliEbyEPidRatioHelper::fgkfHistRangeRap[0],  
170                            AliEbyEPidRatioHelper::fgkfHistRangePhi[0], 
171                            AliEbyEPidRatioHelper::fgkfHistRangePt[0],   
172                            AliEbyEPidRatioHelper::fgkfHistRangeSign[0], 
173                            0.5, -0.5, -3.,-0.5};
174   
175   Double_t maxHnDCA[11] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1],3.5,
176                            AliEbyEPidRatioHelper::fgkfHistRangeEta[1], 
177                            AliEbyEPidRatioHelper::fgkfHistRangeRap[1],  
178                            AliEbyEPidRatioHelper::fgkfHistRangePhi[1], 
179                            AliEbyEPidRatioHelper::fgkfHistRangePt[1],   
180                            AliEbyEPidRatioHelper::fgkfHistRangeSign[1], 
181                            4.5, 1.5, 3., 3.5};
182
183
184   fHnDCA = new THnSparseD("hnDCA", "cent:pid:etaRec:yRec:phiRec:ptRec:sign:contPart:contSign:DCArAccepted:DCAr", 11, binHnDCA, minHnDCA, maxHnDCA);
185   
186   fHnDCA->Sumw2();
187
188   fHnDCA->GetAxis(0)->SetTitle("centrality");                   //  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
189   fHnDCA->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}");   //  0 | 1 | 2 | 3 
190   fHnDCA->GetAxis(2)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
191   fHnDCA->GetAxis(3)->SetTitle("#it{y}_{Rec}");                 //  rapidity  [-0.5, 0.5]
192   fHnDCA->GetAxis(4)->SetTitle("#varphi_{Rec} (rad)");          //  phi  [ 0. , 2Pi]
193   fHnDCA->GetAxis(5)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})");  //  pT   [ 0.2, 2.6]
194   fHnDCA->GetAxis(6)->SetTitle("sign");                         //  -1 | 0 | +1 
195   fHnDCA->GetAxis(7)->SetTitle("contPart");                     //  1  primary | 2 missId | 3 from WeakDecay | 4 p from Material
196   fHnDCA->GetAxis(8)->SetTitle("DCArAccepted");                 //  0 not accepted | 1 accepted 
197   fHnDCA->GetAxis(9)->SetTitle("DCAr");                         //  DCAr [-3, 3]
198
199
200   fHelper->BinLogAxis(fHnDCA,  5, fESDTrackCuts);
201   fHelper->BinLogAxis(fHnDCA,  5, fESDTrackCutsBkg);
202
203   // -- Set binning for DCAr
204   Double_t binsDCAr[77] = {-3.,-2.85,-2.7,-2.55,-2.4,-2.25,-2.1,-1.95,-1.8,-1.65,-1.5,-1.35,-1.2,-1.05,-0.9,-0.75,-0.6,-0.45,-0.3,-0.285,-0.27,-0.255,-0.24,-0.225,-0.21,-0.195,-0.18,-0.165,-0.15,-0.135,-0.12,-0.105,-0.09,-0.075,-0.06,-0.045,-0.03,-0.015,0.,0.015,0.03,0.045,0.06,0.075,0.09,0.105,0.12,0.135,0.15,0.165,0.18,0.195,0.21,0.225,0.24,0.255,0.27,0.285,0.3,0.45,0.6,0.75,0.9,1.05,1.2,1.35,1.5,1.65,1.8,1.95,2.1,2.25,2.4,2.55,2.7,2.85,3.};
205   fHnDCA->GetAxis(9)->Set(76, binsDCAr);
206
207   // ------------------------------------------------------------------
208   
209   return;
210 }
211
212 //________________________________________________________________________
213 Int_t AliEbyEPidRatioDCA::GetContIdxTrack(Int_t label, Int_t sign, Int_t gPdgCode) {
214   Int_t contIdx = -1;
215
216   AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
217   if (!particle)
218     return contIdx;
219
220   Bool_t isPhysicalPrimary = (fESD) ? fStack->IsPhysicalPrimary(label): (static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary();
221   Bool_t isSecondaryFromWeakDecay = (fESD) ? fStack->IsSecondaryFromWeakDecay(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromWeakDecay();
222   Bool_t isSecondaryFromMaterial  = (fESD) ? fStack->IsSecondaryFromMaterial(label)  : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromMaterial();
223   
224   if (isPhysicalPrimary) {
225     if (gPdgCode == 0) {
226       // -- Check if correctly identified 
227       if (particle->PdgCode() == (sign*gPdgCode))
228         contIdx = 1;   
229       // -- MissIdentification
230       else 
231         contIdx = 2;
232     }
233     else
234       contIdx = 1;   
235   }
236
237   // -- Check if secondaries from material or weak decay
238   else if(isSecondaryFromWeakDecay)
239     contIdx = 3;
240   else if (isSecondaryFromMaterial)
241     contIdx = 4;
242   else
243     contIdx = -1;
244
245   return contIdx;
246 }
247
248
249