]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidTTaskMC.cxx
Update PR task: drathee
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidTTaskMC.cxx
CommitLineData
027ebdf2 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"
9e4f5bb3 43#include "AliPIDCombined.h"
027ebdf2 44#include "AliHelperPID.h"
45using std::endl;
46using std::cout;
47#include "AliEbyEPidTTaskMC.h"
48
49ClassImp(AliEbyEPidTTaskMC)
50
51//-----------------------------------------------------------------------
52AliEbyEPidTTaskMC::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),
027ebdf2 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 */
50c07b57 80 DefineInput(0,TChain::Class());
81 DefineOutput(1, TList::Class());
82 DefineOutput(2, TTree::Class());
027ebdf2 83}
84
85AliEbyEPidTTaskMC::~AliEbyEPidTTaskMC() {
86 //! Cleaning up
50c07b57 87 if (fThnList) delete fThnList;
88 if (fHelperPID) delete fHelperPID;
89 if (fEventTree) delete fEventTree;
027ebdf2 90}
91
92//---------------------------------------------------------------------------------
93void 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 }
50c07b57 104 TDirectory *owd = gDirectory;
105 OpenFile(1);
027ebdf2 106 fEventTree = new TTree("fEventTree","fEventTree");
50c07b57 107 owd->cd();
108
027ebdf2 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//----------------------------------------------------------------------------------
137void AliEbyEPidTTaskMC::UserExec( Option_t * ){
138
50c07b57 139 fEventCounter->Fill(1);
140
027ebdf2 141 AliAODEvent* event = dynamic_cast<AliAODEvent*>(InputEvent());
142 if (!event) {
143 Printf("ERROR 01: AOD not found ");
144 return;
145 }
50c07b57 146
027ebdf2 147 fEventCounter->Fill(2);
148
149 Int_t gCent = -1;
638ae862 150 //Float_t gRefMul = -1;
027ebdf2 151
152 AliAODHeader *aodHeader = event->GetHeader();
153 gCent = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
638ae862 154 // gRefMul = aodHeader->GetRefMultiplicity();
027ebdf2 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
122cbc27 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
027ebdf2 192 fEventCounter->Fill(6);
50c07b57 193
027ebdf2 194 fRunNumber = event->GetRunNumber();
195 fFilterBit = fAODtrackCutBit;
196 fCentPercentile = gCent;
197 fVertexX = vertex->GetX();
198 fVertexY = vertex->GetY();
199 fVertexZ = vertex->GetZ();
9e4f5bb3 200
201//Default TPC priors
202 if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
203
027ebdf2 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);
9e4f5bb3 213 // if(a < 0 || a > 2) continue;
027ebdf2 214 fEventCounter->Fill(13);
122cbc27 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;
027ebdf2 221
122cbc27 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
9e4f5bb3 234 // cout << b << " " << track->Charge() << " " << isph << endl;
027ebdf2 235
236 fTrackPt[iTracks] = (Float_t)track->Pt();
237 fTrackPhi[iTracks] = (Float_t)track->Phi();
238 fTrackEta[iTracks] = (Float_t)track->Eta();
122cbc27 239 fTrackCharge[iTracks] = isph;
240 fTrackPid[iTracks] = b;
027ebdf2 241 iTracks++;
242 }
243 fNumberOfTracks = iTracks;
244
122cbc27 245
50c07b57 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);
122cbc27 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
50c07b57 273 //Int_t a = 0;
274 // Int_t a = fHelperPID->GetMCParticleSpecie((AliVEvent*) event,(AliVTrack*)partMC,1);
027ebdf2 275
50c07b57 276 //if(a < 0 || a > 2) continue;
122cbc27 277 // Int_t icharge = a > 0 ? 0 : 1;
278 // cout << a << " " << b << " " << iMC << " " << mTracks << " " << k << endl;
027ebdf2 279
50c07b57 280 fTrackPtM[mTracks] = (Float_t)partMC->Pt();
281 fTrackPhiM[mTracks] = (Float_t)partMC->Phi();
282 fTrackEtaM[mTracks] = (Float_t)partMC->Eta();
122cbc27 283 fTrackChargeM[mTracks] = k;
50c07b57 284 fTrackPidM[mTracks] = a;
285 mTracks++;
286 }
287 // cout << mTracks << " " << nMC << endl;
288 fEventCounter->Fill(26);
289 fNumberOfTracksM = mTracks;
027ebdf2 290
291 fEventCounter->Fill(30);
292 fEventTree->Fill();
293
294 if(isMC) fEventCounter->Fill(46);
295 else fEventCounter->Fill(48);
50c07b57 296
027ebdf2 297 PostData(1, fThnList);
298 PostData(2, fEventTree);
299}
300
301//___________________________________________________________
302void AliEbyEPidTTaskMC::Terminate( Option_t * ){
303 Info("AliEbyEPidTTaskMC"," Task Successfully finished");
304}
305
306//___________________________________________________________
307Bool_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;
122cbc27 311 if (TMath::Abs(track->Eta()) > 0.8) return kFALSE;
027ebdf2 312 return kTRUE;
313}
314
315
316//___________________________________________________________
317Bool_t AliEbyEPidTTaskMC::AcceptMCTrack(AliAODMCParticle *track) const {
318 if(!track) return kFALSE;
027ebdf2 319 if (track->Charge() == 0 ) return kFALSE;
122cbc27 320 if (TMath::Abs(track->Eta()) > 0.8) return kFALSE;
027ebdf2 321 return kTRUE;
322}