78461a8c77b3e2628ea3c03dec8f6128406b6c12
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidTTaskMC.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 //          Dealing with Wide pT Addition (Test Only - not for use)
22 //=========================================================================//
23
24 #include "TChain.h"
25 #include "TList.h"
26 #include "TFile.h"
27 #include "TTree.h"
28 #include "TH1D.h"
29 #include "TH2F.h"
30 #include "TH3F.h"
31 #include "TCanvas.h"
32 #include "AliAnalysisTask.h"
33 #include "AliAnalysisManager.h"
34 #include "AliVEvent.h"
35 #include "AliESD.h"
36 #include "AliESDEvent.h"
37 #include "AliAODEvent.h"
38 #include "AliAODMCParticle.h"
39 #include "AliAODMCHeader.h"
40 #include "AliPIDResponse.h"
41 #include "AliAODHeader.h"
42 #include "AliAODpidUtil.h"
43 #include "AliPIDCombined.h"
44 #include "AliHelperPID.h"
45 using std::endl;
46 using std::cout;
47 #include "AliEbyEPidTTaskMC.h"
48
49 ClassImp(AliEbyEPidTTaskMC)
50
51 //-----------------------------------------------------------------------
52 AliEbyEPidTTaskMC::AliEbyEPidTTaskMC( const char *name ) : AliAnalysisTaskSE( name ), 
53   fThnList(NULL), 
54   isMC(kFALSE), 
55   fCentralityEstimator("V0M"), 
56   fAODtrackCutBit(128),
57   fHelperPID(0x0),
58   fEventCounter(NULL), 
59   fEventTree(NULL),
60   fRunNumber(0),
61   fNumberOfTracks(8000),
62   fNumberOfTracksM(8000),
63   fFilterBit(99),
64   fCentPercentile(-1),
65   fVertexX(-99),
66   fVertexY(-99),
67   fVertexZ(-99) {
68   /*
69     fTrackPt[kTrack];
70     fTrackEta[kTrack];
71     fTrackPhi[kTrack];
72     fTrackPtM[kTrack];
73     fTrackEtaM[kTrack];
74     fTrackPhiM[kTrack];
75     fTrackCharge[kTrack];
76     fTrackPid[kTrack];
77     fTrackChargeM[kTrack];
78     fTrackPidM[kTrack];
79   */
80   DefineInput(0,TChain::Class());
81   DefineOutput(1, TList::Class()); 
82   DefineOutput(2, TTree::Class()); 
83 }
84
85 AliEbyEPidTTaskMC::~AliEbyEPidTTaskMC() {
86   //!   Cleaning up
87    if (fThnList)   delete fThnList;
88    if (fHelperPID) delete fHelperPID;
89    if (fEventTree) delete fEventTree;
90 }
91
92 //---------------------------------------------------------------------------------
93 void AliEbyEPidTTaskMC::UserCreateOutputObjects() {
94   fThnList = new TList();
95   fThnList->SetOwner(kTRUE);
96
97   fEventCounter = new TH1D("fEventCounter","EventCounter", 70, 0.5,70.5);
98   fThnList->Add(fEventCounter);
99
100   TList *ll = (TList*)fHelperPID->GetOutputList();
101   for (Int_t ikey = 0; ikey < ll->GetEntries(); ikey++) {
102     fThnList->Add(ll->At(ikey));
103   }
104   TDirectory *owd = gDirectory;
105   OpenFile(1);  
106   fEventTree = new TTree("fEventTree","fEventTree");
107   owd->cd();
108  
109   fEventTree->Branch("fRunNumber",      &fRunNumber,     "fRunNumber/I");
110   fEventTree->Branch("fFilterBit",      &fFilterBit,     "fFilterBit/I");
111   fEventTree->Branch("fNumberOfTracks", &fNumberOfTracks,"fNumberOfTracks/I");
112   fEventTree->Branch("fCentPercentile", &fCentPercentile,"fCentPercentile/F");
113   
114   fEventTree->Branch("fVertexX",        &fVertexX,       "fVertexX/F");
115   fEventTree->Branch("fVertexY",        &fVertexY,       "fVertexY/F");
116   fEventTree->Branch("fVertexZ",        &fVertexZ,       "fVertexZ/F");
117   
118   fEventTree->Branch("fTrackPt",        fTrackPt,        "fTrackPt[fNumberOfTracks]/F");
119   fEventTree->Branch("fTrackPhi",       fTrackPhi,       "fTrackPhi[fNumberOfTracks]/F");
120   fEventTree->Branch("fTrackEta",       fTrackEta,       "fTrackEta[fNumberOfTracks]/F");
121   fEventTree->Branch("fTrackCharge",    fTrackCharge,    "fTrackCharge[fNumberOfTracks]/I");
122   fEventTree->Branch("fTrackPid",       fTrackPid,       "fTrackPid[fNumberOfTracks]/I");
123   
124
125   fEventTree->Branch("fNumberOfTracksM", &fNumberOfTracksM, "fNumberOfTracksM/I");
126   fEventTree->Branch("fTrackPtM",        fTrackPtM,        "fTrackPtM[fNumberOfTracksM]/F");
127   fEventTree->Branch("fTrackPhiM",       fTrackPhiM,       "fTrackPhiM[fNumberOfTracksM]/F");
128   fEventTree->Branch("fTrackEtaM",       fTrackEtaM,       "fTrackEtaM[fNumberOfTracksM]/F");
129   fEventTree->Branch("fTrackChargeM",    fTrackChargeM,    "fTrackChargeM[fNumberOfTracksM]/I");
130   fEventTree->Branch("fTrackPidM",       fTrackPidM,       "fTrackPidM[fNumberOfTracksM]/I");
131  
132   PostData(1, fThnList);
133   PostData(2, fEventTree);  
134 }
135
136 //----------------------------------------------------------------------------------
137 void AliEbyEPidTTaskMC::UserExec( Option_t * ){
138
139   fEventCounter->Fill(1);
140    
141   AliAODEvent* event = dynamic_cast<AliAODEvent*>(InputEvent());
142   if (!event) {
143     Printf("ERROR 01: AOD not found ");
144     return;
145   }
146   
147   fEventCounter->Fill(2);
148
149   Int_t gCent   = -1;
150   //Float_t gRefMul = -1;
151   
152   AliAODHeader *aodHeader = event->GetHeader();
153   gCent = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
154  // gRefMul = aodHeader->GetRefMultiplicity();
155   if (gCent < 0 || gCent > 100) return;
156   fEventCounter->Fill(3);  
157
158   const AliAODVertex *vertex = event->GetPrimaryVertex();
159   if(!vertex) return;
160   fEventCounter->Fill(4);
161   Bool_t vtest = kFALSE;
162   Double32_t fCov[6];
163   vertex->GetCovarianceMatrix(fCov);
164   if(vertex->GetNContributors() > 0) {
165       if(fCov[5] != 0) {
166         vtest = kTRUE;
167       }
168   }
169   if(!vtest)return;
170   
171   fEventCounter->Fill(5);
172
173   AliCentrality *centrality = event->GetCentrality();
174   if (centrality->GetQuality() != 0) return;
175
176  fEventCounter->Fill(21);
177   
178   TClonesArray *arrayMC= 0; 
179   arrayMC = dynamic_cast<TClonesArray*> (event->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
180   if (!arrayMC) {
181     Printf("Error: MC particles branch not found!\n");
182       return;
183   }
184   fEventCounter->Fill(22);
185   AliAODMCHeader *mcHdr=0;
186   mcHdr=(AliAODMCHeader*)event->GetList()->FindObject(AliAODMCHeader::StdBranchName());  
187   if(!mcHdr) {
188       Printf("MC header branch not found!\n");
189       return;
190   }
191
192   fEventCounter->Fill(6);
193   
194   fRunNumber = event->GetRunNumber();
195   fFilterBit = fAODtrackCutBit;
196   fCentPercentile = gCent;
197   fVertexX = vertex->GetX();
198   fVertexY = vertex->GetY();
199   fVertexZ = vertex->GetZ();
200  
201 //Default TPC priors
202   if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
203  
204   Int_t iTracks = 0; 
205   for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
206     AliAODTrack* track = dynamic_cast<AliAODTrack *>(event->GetTrack(itrk));
207     fEventCounter->Fill(10);
208     if (!track) continue;
209     fEventCounter->Fill(11);
210     if (!AcceptTrack(track)) continue;
211     fEventCounter->Fill(12);
212     Int_t a = fHelperPID->GetParticleSpecies((AliVTrack*) track,kTRUE);
213     // if(a < 0 || a > 2) continue;
214     fEventCounter->Fill(13);
215     //Int_t icharge = track->Charge() > 0 ? 0 : 1;
216     Int_t b = -999;
217     if (a == 0 ) b = 1;
218     else if (a == 1 ) b = 2;
219     else if (a == 2 ) b = 3;
220     else b = 4;
221
222     if (track->Charge()  < 0 ) b = -1*b;
223    
224     Int_t isph=-999;
225     if (arrayMC) {
226       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
227       if (!partMC) {
228         AliError("Cannot get MC particle");
229         continue;
230       }
231       isph=partMC->IsPhysicalPrimary();
232     }
233
234     //   cout << b << "  " << track->Charge() << "  " << isph << endl;
235
236     fTrackPt[iTracks]     = (Float_t)track->Pt();
237     fTrackPhi[iTracks]    = (Float_t)track->Phi();
238     fTrackEta[iTracks]    = (Float_t)track->Eta();
239     fTrackCharge[iTracks] = isph;
240     fTrackPid[iTracks] = b;
241     iTracks++;
242   }
243   fNumberOfTracks = iTracks;
244   
245    
246   fEventCounter->Fill(23);
247   
248   
249   Int_t mTracks = 0; 
250
251   Int_t nMC = arrayMC->GetEntries();
252   
253   for (Int_t iMC = 0; iMC < nMC; iMC++) {
254     //fEventCounter->Fill(24);
255     AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
256     // if(!partMC) continue;
257     if(!AcceptMCTrack(partMC)) continue;
258     // fEventCounter->Fill(25);
259     Int_t b  = partMC->GetPdgCode();
260     Int_t k = -99;
261     k = partMC->IsPhysicalPrimary();
262    
263     Int_t a = -999;
264     if (b == 211 )   a =  1;
265     else if (b == -211 )  a = -1;
266     else if (b == 2212 )  a =  3;
267     else if (b == -2212 ) a = -3;
268     else if (b == 321 )   a =  2;
269     else if (b == -321 )  a = -2;
270     else if (b > 0 )  a = 4;
271     else if (b < 0 )  a = -4;
272
273     //Int_t a = 0;      
274     // Int_t a = fHelperPID->GetMCParticleSpecie((AliVEvent*) event,(AliVTrack*)partMC,1);
275     
276     //if(a < 0 || a > 2) continue;
277     // Int_t icharge = a > 0 ? 0 : 1;
278     //  cout << a << "  " << b << "  " << iMC << "  " << mTracks << "  " << k << endl;
279     
280     fTrackPtM[mTracks]     = (Float_t)partMC->Pt();
281     fTrackPhiM[mTracks]    = (Float_t)partMC->Phi();
282     fTrackEtaM[mTracks]    = (Float_t)partMC->Eta();
283     fTrackChargeM[mTracks] = k;
284     fTrackPidM[mTracks] = a;
285     mTracks++;
286   }
287   //  cout << mTracks << "  " << nMC << endl;
288   fEventCounter->Fill(26);
289   fNumberOfTracksM = mTracks;
290   
291   fEventCounter->Fill(30);
292   fEventTree->Fill();
293   
294   if(isMC) fEventCounter->Fill(46);
295   else fEventCounter->Fill(48);
296      
297   PostData(1, fThnList);  
298   PostData(2, fEventTree);
299 }
300
301 //___________________________________________________________
302 void AliEbyEPidTTaskMC::Terminate( Option_t * ){
303   Info("AliEbyEPidTTaskMC"," Task Successfully finished");
304 }
305
306 //___________________________________________________________
307 Bool_t AliEbyEPidTTaskMC::AcceptTrack(AliAODTrack *track) const {
308   if(!track)                                  return kFALSE;
309   if (track->Charge() == 0 )                  return kFALSE;
310   if (!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;
311   if (TMath::Abs(track->Eta()) > 0.8) return kFALSE;
312   return kTRUE;
313 }
314
315
316 //___________________________________________________________
317 Bool_t AliEbyEPidTTaskMC::AcceptMCTrack(AliAODMCParticle *track) const {
318   if(!track) return kFALSE;
319   if (track->Charge() == 0 ) return kFALSE;
320   if (TMath::Abs(track->Eta()) > 0.8) return kFALSE;
321   return kTRUE;
322 }