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