]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidTTask.cxx
AliAODEvent::GetHeader now return AliVHeader
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidTTask.cxx
CommitLineData
3435d75f 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"
e69ef551 41#include "AliPIDCombined.h"
3435d75f 42#include "AliAODHeader.h"
43#include "AliAODpidUtil.h"
44#include "AliHelperPID.h"
45using std::endl;
46using std::cout;
47#include "AliEbyEPidTTask.h"
48
49ClassImp(AliEbyEPidTTask)
50
51//-----------------------------------------------------------------------
52AliEbyEPidTTask::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
87AliEbyEPidTTask::~AliEbyEPidTTask() {
88 //! Cleaning up
89 if (fThnList) delete fThnList;
90 if (fHelperPID) delete fHelperPID;
e69ef551 91 if (fEventTree) delete fEventTree;
3435d75f 92}
93
94//---------------------------------------------------------------------------------
95void AliEbyEPidTTask::UserCreateOutputObjects() {
96 fThnList = new TList();
97 fThnList->SetOwner(kTRUE);
98
be2d8a61 99 fEventCounter = new TH1D("fEventCounter","EventCounter", 70, 0.5,70.5);
3435d75f 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 }
e69ef551 106 TDirectory *owd = gDirectory;
107 OpenFile(1);
3435d75f 108 fEventTree = new TTree("fEventTree","fEventTree");
e69ef551 109 owd->cd();
110
3435d75f 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//----------------------------------------------------------------------------------
140void 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
be2d8a61 150 fEventCounter->Fill(2);
151
3435d75f 152 Int_t gCent = -1;
638ae862 153 // Float_t gRefMul = -1;
3435d75f 154
155 AliAODHeader *aodHeader = event->GetHeader();
156 gCent = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
638ae862 157 // gRefMul = aodHeader->GetRefMultiplicity();
3435d75f 158 if (gCent < 0 || gCent > 100) return;
be2d8a61 159 fEventCounter->Fill(3);
3435d75f 160
161 const AliAODVertex *vertex = event->GetPrimaryVertex();
162 if(!vertex) return;
be2d8a61 163 fEventCounter->Fill(4);
3435d75f 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
be2d8a61 174 fEventCounter->Fill(5);
175
3435d75f 176 AliCentrality *centrality = event->GetCentrality();
177 if (centrality->GetQuality() != 0) return;
178
be2d8a61 179 fEventCounter->Fill(6);
180
3435d75f 181 fRunNumber = event->GetRunNumber();
182 fFilterBit = fAODtrackCutBit;
183 fCentPercentile = gCent;
184 fVertexX = vertex->GetX();
185 fVertexY = vertex->GetY();
186 fVertexZ = vertex->GetZ();
187
e69ef551 188//Default TPC priors
189 if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
190
191
192
193
3435d75f 194 Int_t iTracks = 0;
195 for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
196 AliAODTrack* track = dynamic_cast<AliAODTrack *>(event->GetTrack(itrk));
be2d8a61 197 fEventCounter->Fill(10);
3435d75f 198 if (!track) continue;
be2d8a61 199 fEventCounter->Fill(11);
3435d75f 200 if (!AcceptTrack(track)) continue;
be2d8a61 201 fEventCounter->Fill(12);
3435d75f 202 Int_t a = fHelperPID->GetParticleSpecies((AliVTrack*) track,kTRUE);
122cbc27 203 // if(a < 0 || a > 2) continue;
be2d8a61 204 fEventCounter->Fill(13);
122cbc27 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
e69ef551 214 // Int_t icharge = track->Charge() > 0 ? 0 : 1;
3435d75f 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();
e69ef551 221 fTrackCharge[iTracks] = track->Charge();
122cbc27 222 fTrackPid[iTracks] = b;
3435d75f 223 iTracks++;
224 }
225 fNumberOfTracks = iTracks;
226
227 if(isMC) {
be2d8a61 228 fEventCounter->Fill(21);
3435d75f 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 }
be2d8a61 235 fEventCounter->Fill(22);
3435d75f 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 }
122cbc27 242
be2d8a61 243 fEventCounter->Fill(23);
3435d75f 244
122cbc27 245 Int_t nMC = arrayMC->GetEntries();
246 Int_t mTracks = 0;
3435d75f 247 for (Int_t iMC = 0; iMC < nMC; iMC++) {
be2d8a61 248 fEventCounter->Fill(24);
122cbc27 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++;
3435d75f 269 }
be2d8a61 270 fEventCounter->Fill(26);
3435d75f 271 fNumberOfTracksM = mTracks;
272 }
273
be2d8a61 274 fEventCounter->Fill(30);
275 fEventTree->Fill();
3435d75f 276
277 PostData(1, fThnList);
278 PostData(2, fEventTree);
279}
280
281//___________________________________________________________
282void AliEbyEPidTTask::Terminate( Option_t * ){
283 Info("AliEbyEPidTTask"," Task Successfully finished");
284}
285
286//___________________________________________________________
287Bool_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//___________________________________________________________
296Bool_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}