]>
Commit | Line | Data |
---|---|---|
980821ba | 1 | // $Id$ |
f52aa4d5 | 2 | // |
3 | // | |
4 | // | |
5 | // | |
6 | ||
33d8da67 | 7 | #include "TChain.h" |
8 | #include "TTree.h" | |
9 | #include "TH1F.h" | |
10 | #include "TH2F.h" | |
11 | #include "TCanvas.h" | |
12 | ||
13 | #include "AliAnalysisTask.h" | |
14 | #include "AliAnalysisManager.h" | |
15 | ||
16 | #include "AliESDEvent.h" | |
17 | #include "AliESDHeader.h" | |
18 | #include "AliESDUtils.h" | |
19 | #include "AliESDInputHandler.h" | |
20 | #include "AliCentrality.h" | |
21 | #include "AliESDpid.h" | |
22 | #include "AliKFParticle.h" | |
23 | ||
24 | #include "AliMCEventHandler.h" | |
25 | #include "AliMCEvent.h" | |
26 | #include "AliStack.h" | |
27 | #include "TParticle.h" | |
7a7744db | 28 | #include "AliAODMCParticle.h" |
33d8da67 | 29 | |
30 | ||
31 | #include "AliESDtrackCuts.h" | |
32 | #include "AliESDv0.h" | |
33 | #include "AliV0vertexer.h" | |
34 | #include "AliESDCaloCluster.h" | |
35 | #include "AliESDCaloCells.h" | |
7a7744db | 36 | #include "AliAODEvent.h" |
33d8da67 | 37 | #include "AliEMCALGeometry.h" |
38 | #include "AliEMCALRecoUtils.h" | |
7a7744db | 39 | #include "AliOADBContainer.h" |
40 | #include "AliVEvent.h" | |
41 | #include "AliVVertex.h" | |
33d8da67 | 42 | #include "AliVCluster.h" |
7a7744db | 43 | #include "AliVCaloCells.h" |
33d8da67 | 44 | #include "AliAnalysisTaskEMCALClusterizeFast.h" |
45 | #include "TLorentzVector.h" | |
46 | ||
47 | ||
48 | #include "AliAnalysisTaskEMCALPhoton.h" | |
49 | #include "TFile.h" | |
50 | ||
51 | ||
c64cb1f6 | 52 | using std::cout; |
53 | using std::endl; | |
54 | ||
33d8da67 | 55 | ClassImp(AliAnalysisTaskEMCALPhoton) |
56 | ||
57 | //________________________________________________________________________ | |
f1c66719 | 58 | AliAnalysisTaskEMCALPhoton::AliAnalysisTaskEMCALPhoton() : |
59 | AliAnalysisTaskSE(), | |
60 | fTrCuts(0), | |
61 | fPrTrCuts(0), | |
62 | fSelTracks(0), | |
63 | fSelPrimTracks(0), | |
7a7744db | 64 | fTracks(0), |
f1c66719 | 65 | fPhotConvArray(0), |
66 | fMyClusts(0), | |
67 | fMyAltClusts(0), | |
68 | fMyCells(0), | |
69 | fMyTracks(0), | |
70 | fMyMcParts(0), | |
71 | fHeader(0x0), | |
5b19b2cc | 72 | fOADBContainer(0), |
f1c66719 | 73 | fCaloClusters(0), |
74 | fCaloClustersNew(0), | |
7a7744db | 75 | fAODMCParticles(0), |
76 | fVCells(0), | |
f1c66719 | 77 | fGeom(0x0), |
78 | fTimeResTOF(0), | |
79 | fMipResponseTPC(0), | |
80 | fGeoName("EMCAL_COMPLETEV1"), | |
81 | fPeriod("LHC11d"), | |
82 | fIsTrain(0), | |
83 | fIsMC(0), | |
40b40be6 | 84 | fDebug(0), |
c5f3236f | 85 | fRedoV0(1), |
f1c66719 | 86 | fIsGrid(0), |
87 | fClusThresh(2.0), | |
88 | fClusterizer(0), | |
89 | fCaloClustersName("EmcalClusterv2"), | |
90 | fESD(0), | |
7a7744db | 91 | fAOD(0), |
92 | fVev(0), | |
f1c66719 | 93 | fMCEvent(0), |
94 | fStack(0x0), | |
95 | fOutputList(0), | |
96 | fTree(0), | |
c5f3236f | 97 | fMyMcIndex(0), |
f1c66719 | 98 | fNV0sBefAndAftRerun(0), |
99 | fConversionVtxXY(0), | |
100 | fInvMassV0(0), | |
101 | fInvMassV0KF(0), | |
102 | fInvMassV0SS(0), | |
103 | fDedxPAll(0) | |
104 | { | |
105 | // Default constructor. | |
7a7744db | 106 | for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0; |
f1c66719 | 107 | } |
108 | ||
109 | //________________________________________________________________________ | |
110 | AliAnalysisTaskEMCALPhoton::AliAnalysisTaskEMCALPhoton(const char *name) : | |
111 | AliAnalysisTaskSE(name), | |
33d8da67 | 112 | fTrCuts(0), |
113 | fPrTrCuts(0), | |
114 | fSelTracks(0), | |
115 | fSelPrimTracks(0), | |
7a7744db | 116 | fTracks(0), |
33d8da67 | 117 | fPhotConvArray(0), |
118 | fMyClusts(0), | |
119 | fMyAltClusts(0), | |
120 | fMyCells(0), | |
121 | fMyTracks(0), | |
122 | fMyMcParts(0), | |
123 | fHeader(0), | |
5b19b2cc | 124 | fOADBContainer(0), |
33d8da67 | 125 | fCaloClusters(0), |
126 | fCaloClustersNew(0), | |
7a7744db | 127 | fAODMCParticles(0), |
128 | fVCells(0), | |
33d8da67 | 129 | fGeom(0x0), |
130 | fTimeResTOF(0), | |
131 | fMipResponseTPC(0), | |
132 | fGeoName("EMCAL_COMPLETEV1"), | |
133 | fPeriod("LHC11c"), | |
134 | fIsTrain(0), | |
135 | fIsMC(0), | |
40b40be6 | 136 | fDebug(0), |
c5f3236f | 137 | fRedoV0(1), |
33d8da67 | 138 | fIsGrid(0), |
139 | fClusThresh(2.), | |
140 | fClusterizer(0), | |
141 | fCaloClustersName("EmcalClusterv2"), | |
33d8da67 | 142 | fESD(0), |
7a7744db | 143 | fAOD(0), |
144 | fVev(0), | |
33d8da67 | 145 | fMCEvent(0), |
146 | fStack(0x0), | |
33d8da67 | 147 | fOutputList(0), |
148 | fTree(0), | |
c5f3236f | 149 | fMyMcIndex(0), |
33d8da67 | 150 | fNV0sBefAndAftRerun(0), |
151 | fConversionVtxXY(0), | |
152 | fInvMassV0(0), | |
153 | fInvMassV0KF(0), | |
154 | fInvMassV0SS(0), | |
155 | fDedxPAll(0) | |
33d8da67 | 156 | { |
157 | // Constructor | |
158 | ||
159 | // Define input and output slots here | |
160 | // Input slot #0 works with a TChain | |
161 | DefineInput(0, TChain::Class()); | |
162 | // Output slot #0 id reserved by the base class for AOD | |
163 | // Output slot #1 writes into a TH1 container | |
164 | DefineOutput(1, TList::Class()); | |
ec681cf3 | 165 | DefineOutput(2, TTree::Class()); |
33d8da67 | 166 | } |
f1c66719 | 167 | |
33d8da67 | 168 | //________________________________________________________________________ |
169 | void AliAnalysisTaskEMCALPhoton::UserCreateOutputObjects() | |
170 | { | |
f1c66719 | 171 | // Create histograms, called once. |
4b6ab175 | 172 | if(this->fDebug) |
c5f3236f | 173 | printf("::UserCreateOutputObjects() starting\n"); |
174 | ||
33d8da67 | 175 | fSelTracks = new TObjArray(); |
176 | ||
177 | fSelPrimTracks = new TObjArray(); | |
178 | ||
b2d49404 | 179 | if (TClass::GetClass("AliPhotonHeaderObj")) |
180 | TClass::GetClass("AliPhotonHeaderObj")->IgnoreTObjectStreamer(); | |
181 | fHeader = new AliPhotonHeaderObj; | |
33d8da67 | 182 | |
b2d49404 | 183 | if (TClass::GetClass("AliPhotonConvObj")) |
184 | TClass::GetClass("AliPhotonConvObj")->IgnoreTObjectStreamer(); | |
185 | fPhotConvArray = new TClonesArray("AliPhotonConvObj"); | |
33d8da67 | 186 | |
b2d49404 | 187 | if (TClass::GetClass("AliPhotonClusterObj")) |
188 | TClass::GetClass("AliPhotonClusterObj")->IgnoreTObjectStreamer(); | |
189 | fMyClusts = new TClonesArray("AliPhotonClusterObj"); | |
33d8da67 | 190 | |
b2d49404 | 191 | if (TClass::GetClass("AliPhotonCellObj")) |
192 | TClass::GetClass("AliPhotonCellObj")->IgnoreTObjectStreamer(); | |
193 | fMyCells = new TClonesArray("AliPhotonCellObj"); | |
33d8da67 | 194 | |
b2d49404 | 195 | if (TClass::GetClass("AliPhotonTrackObj")) |
196 | TClass::GetClass("AliPhotonTrackObj")->IgnoreTObjectStreamer(); | |
197 | fMyTracks = new TClonesArray("AliPhotonTrackObj"); | |
33d8da67 | 198 | |
199 | if (fClusterizer || fIsTrain){ | |
200 | if(fClusterizer) | |
201 | fCaloClustersName = fClusterizer->GetNewClusterArrayName(); | |
202 | else { | |
203 | if(fPeriod.Contains("10h") || fPeriod.Contains("11h")) | |
204 | fCaloClustersName = "EmcalClustersv1"; | |
205 | else | |
206 | fCaloClustersName = "EmcalClustersv2"; | |
207 | } | |
b2d49404 | 208 | if (TClass::GetClass("AliPhotonClusterObj")) |
209 | TClass::GetClass("AliPhotonClusterObj")->IgnoreTObjectStreamer(); | |
210 | fMyAltClusts = new TClonesArray("AliPhotonClusterObj"); | |
33d8da67 | 211 | } |
212 | cout<<fCaloClustersName<<endl; | |
213 | if(fIsMC){ | |
b2d49404 | 214 | if (TClass::GetClass("AliPhotonMcPartObj")) |
215 | TClass::GetClass("AliPhotonMcPartObj")->IgnoreTObjectStreamer(); | |
216 | fMyMcParts = new TClonesArray("AliPhotonMcPartObj"); | |
33d8da67 | 217 | } |
218 | ||
7a7744db | 219 | fCaloClusters = new TClonesArray(); |
33d8da67 | 220 | |
221 | fOutputList = new TList(); | |
222 | fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging) | |
223 | ||
224 | if( !fTree){ | |
55dcd0ae | 225 | TFile *f = OpenFile(2); |
33d8da67 | 226 | TDirectory::TContext context(f); |
227 | if (f) { | |
228 | f->SetCompressionLevel(2); | |
229 | fTree = new TTree("photon_ana_out", "StandaloneTree"); | |
230 | fTree->SetDirectory(f); | |
231 | if (fIsTrain) { | |
232 | fTree->SetAutoFlush(-2*1024*1024); | |
233 | fTree->SetAutoSave(0); | |
234 | } else { | |
235 | fTree->SetAutoFlush(-32*1024*1024); | |
236 | fTree->SetAutoSave(0); | |
237 | } | |
238 | ||
239 | fTree->Branch("header", &fHeader, 16*1024, 99); | |
240 | fTree->Branch("conversions", &fPhotConvArray, 8*16*1024, 99); | |
241 | fTree->Branch("clusters", &fMyClusts, 8*16*1024, 99); | |
242 | if(fClusterizer || fIsTrain) | |
243 | fTree->Branch(fCaloClustersName, &fMyAltClusts, 8*16*1024, 99); | |
244 | fTree->Branch("cells", &fMyCells, 8*16*1024, 99); | |
3b37c011 | 245 | fTree->Branch("IsoTracks", &fMyTracks, 8*16*1024, 99); |
33d8da67 | 246 | if(fIsMC) |
247 | fTree->Branch("mcparts", &fMyMcParts, 8*16*1024, 99); | |
248 | } | |
249 | } | |
55dcd0ae | 250 | //if(fIsGrid)fOutputList->Add(fTree); |
33d8da67 | 251 | fGeom = AliEMCALGeometry::GetInstance(fGeoName); |
5b19b2cc | 252 | fOADBContainer = new AliOADBContainer("AliEMCALgeo"); |
253 | fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo"); | |
33d8da67 | 254 | |
255 | ||
256 | fNV0sBefAndAftRerun = new TH2F("hNV0sBefAndAftRerun","check if the number of v0s change with rerun;old v0 n;new v0 n",50,0.5,50.5,50,0.5,50.5); | |
257 | fOutputList->Add(fNV0sBefAndAftRerun); | |
258 | ||
259 | fConversionVtxXY = new TH2F("hConversionVtxXY","x and y of conversion vertex candidates;x;y",1000,-100,100,1000,-100,100); | |
260 | fOutputList->Add(fConversionVtxXY); | |
261 | ||
262 | fInvMassV0 = new TH1F("hInvMassV0","v0->GetEffMass();v0->GetEffMass();dN/dM",400,0,4); | |
263 | fOutputList->Add(fInvMassV0); | |
264 | ||
265 | fInvMassV0KF = new TH1F("hInvMassV0KF","Inv. mass calculated from AliKFParticle made of V0 tracks;mass_{TrTr};dN/dM",400,0,4); | |
266 | fOutputList->Add(fInvMassV0KF); | |
267 | ||
268 | fInvMassV0SS = new TH1F("hInvMassV0SS","Inv. mass (same sign) calculated from AliKFParticle made of V0 tracks;mass_{TrTr};dN/dM",400,0,4); | |
269 | fOutputList->Add(fInvMassV0SS); | |
270 | ||
271 | fDedxPAll = new TH2F("hDedxPAll","dE/dx vs p (all selected tracks);p (GeV/c);dE/dx (a.u.)",400,0,40, 150,0,150); | |
272 | fOutputList->Add(fDedxPAll); | |
273 | ||
274 | ||
275 | PostData(1, fOutputList); | |
af1750ae | 276 | PostData(2, fTree); |
c5f3236f | 277 | |
4b6ab175 | 278 | if(this->fDebug) |
c5f3236f | 279 | printf("::UserCreateOutputObjects() DONE!\n"); |
280 | ||
33d8da67 | 281 | } |
282 | ||
283 | //________________________________________________________________________ | |
284 | void AliAnalysisTaskEMCALPhoton::UserExec(Option_t *) | |
285 | { | |
f1c66719 | 286 | // User exec, called once per event. |
287 | ||
d1faab8e | 288 | Bool_t isSelected = kTRUE; |
33d8da67 | 289 | if(fPeriod.Contains("11")){ |
290 | if(fPeriod.Contains("11a")) | |
291 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1); | |
292 | if(fPeriod.Contains("11c") ||fPeriod.Contains("11d") ) | |
293 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7); | |
294 | if(fPeriod.Contains("11h") ) | |
772f286b | 295 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);//kEMCEGA); |
33d8da67 | 296 | |
297 | } | |
c5ac19f6 | 298 | if(fIsMC){ |
299 | isSelected = kTRUE; | |
4b6ab175 | 300 | if(this->fDebug) |
7a7744db | 301 | printf("+++Message+++: MC input files.\n"); |
c5ac19f6 | 302 | } |
7a7744db | 303 | if(!isSelected){ |
4b6ab175 | 304 | if(this->fDebug) |
7a7744db | 305 | printf("+++Message+++: Event didn't pass the selection\n"); |
4a4936e3 | 306 | return; |
7a7744db | 307 | } |
4b6ab175 | 308 | if(this->fDebug){ |
c5f3236f | 309 | printf("::UserExec(): event accepted\n"); |
310 | if(fIsMC) | |
311 | printf("\t in MC mode\n"); | |
312 | } | |
7a7744db | 313 | |
314 | TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree(); | |
315 | TFile *inpfile = (TFile*)tree->GetCurrentFile(); | |
316 | ||
317 | fVev = (AliVEvent*)InputEvent(); | |
318 | if (!fVev) { | |
319 | printf("ERROR: event not available\n"); | |
33d8da67 | 320 | return; |
321 | } | |
7a7744db | 322 | Int_t runnumber = InputEvent()->GetRunNumber() ; |
4b6ab175 | 323 | if(this->fDebug) |
7a7744db | 324 | printf("run number = %d\n",runnumber); |
325 | fESD = dynamic_cast<AliESDEvent*>(fVev); | |
326 | if(!fESD){ | |
327 | fAOD = dynamic_cast<AliAODEvent*>(fVev); | |
328 | if(!fAOD){ | |
329 | printf("ERROR: Invalid type of event!!!\n"); | |
330 | return; | |
331 | } | |
4b6ab175 | 332 | else if(this->fDebug) |
7a7744db | 333 | printf("AOD EVENT!\n"); |
334 | } | |
33d8da67 | 335 | |
7a7744db | 336 | AliVVertex *pv = (AliVVertex*)fVev->GetPrimaryVertex(); |
337 | Bool_t pvStatus = kTRUE; | |
338 | if(fESD){ | |
339 | AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex(); | |
340 | pvStatus = esdv->GetStatus(); | |
341 | } | |
342 | ||
343 | if(!pv) { | |
344 | printf("Error: no primary vertex found!\n"); | |
33d8da67 | 345 | return; |
7a7744db | 346 | } |
33d8da67 | 347 | if(TMath::Abs(pv->GetZ())>15) |
348 | return; | |
4b6ab175 | 349 | if(this->fDebug) |
7a7744db | 350 | printf("+++Message+++: event passed the vertex cut\n"); |
33d8da67 | 351 | |
7a7744db | 352 | if(fESD) |
353 | fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks")); | |
354 | if(fAOD) | |
355 | fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks()); | |
c5ac19f6 | 356 | |
7a7744db | 357 | if(!fTracks){ |
358 | AliError("Track array in event is NULL!"); | |
4b6ab175 | 359 | if(this->fDebug) |
7a7744db | 360 | printf("returning due to not finding tracks!\n"); |
361 | return; | |
362 | } | |
363 | ||
364 | const Int_t Ntracks = fTracks->GetEntriesFast(); | |
33d8da67 | 365 | // Track loop to fill a pT spectrum |
7a7744db | 366 | for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) { |
367 | AliVTrack *track = (AliVTrack*)fTracks->At(iTracks); | |
33d8da67 | 368 | if (!track) |
369 | continue; | |
7a7744db | 370 | AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track); |
371 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track); | |
372 | if (esdTrack){ | |
373 | if (fTrCuts && fTrCuts->IsSelected(track)){ | |
374 | fSelTracks->Add(track); | |
375 | fDedxPAll->Fill(track->P(), track->GetTPCsignal()); | |
376 | } | |
377 | if (fPrTrCuts && fPrTrCuts->IsSelected(track)) | |
378 | fSelPrimTracks->Add(track); | |
33d8da67 | 379 | } |
7a7744db | 380 | else if(aodTrack) |
33d8da67 | 381 | fSelPrimTracks->Add(track); |
7a7744db | 382 | }//track loop |
383 | ||
384 | ||
f1c66719 | 385 | |
c5ac19f6 | 386 | fHeader->fInputFileName = inpfile->GetName(); |
ea7d3479 | 387 | fHeader->fRunNumber = runnumber; |
7a7744db | 388 | fHeader->fTrClassMask = fVev->GetHeader()->GetTriggerMask(); |
389 | fHeader->fTrCluster = fVev->GetHeader()->GetTriggerCluster(); | |
33d8da67 | 390 | AliCentrality *cent = InputEvent()->GetCentrality(); |
391 | fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M"); | |
392 | fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1"); | |
393 | fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK"); | |
7a7744db | 394 | |
5b19b2cc | 395 | TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices"); |
7a7744db | 396 | for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ |
397 | if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3) | |
398 | break; | |
399 | /*if(fESD) | |
33d8da67 | 400 | fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod); |
7a7744db | 401 | else*/ |
402 | // if(event->GetEMCALMatrix(mod)) | |
403 | fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod); | |
404 | fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod); | |
405 | } | |
406 | Int_t trackMult = 0; | |
407 | if(fESD){ | |
408 | AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts(); | |
409 | trackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5); | |
4b6ab175 | 410 | if(this->fDebug) |
7a7744db | 411 | printf("ESD Track mult= %d\n",trackMult); |
412 | } | |
413 | else if(fAOD){ | |
414 | trackMult = Ntracks; | |
4b6ab175 | 415 | if(this->fDebug) |
7a7744db | 416 | printf("AOD Track mult= %d\n",trackMult); |
417 | } | |
418 | if(fESD){ | |
419 | TList *l = fESD->GetList(); | |
420 | fCaloClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters")); | |
421 | fVCells = fESD->GetEMCALCells(); | |
422 | fHeader->fNCells = fVCells->GetNumberOfCells(); | |
4b6ab175 | 423 | if(this->fDebug) |
7a7744db | 424 | printf("ESD cluster mult= %d\n",fCaloClusters->GetEntriesFast()); |
425 | } | |
426 | else if(fAOD){ | |
427 | fCaloClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters()); | |
428 | fVCells = fAOD->GetEMCALCells(); | |
429 | fHeader->fNCells = fVCells->GetNumberOfCells(); | |
4b6ab175 | 430 | if(this->fDebug) |
7a7744db | 431 | printf("AOD cluster mult= %d\n",fCaloClusters->GetEntriesFast()); |
33d8da67 | 432 | } |
7a7744db | 433 | |
434 | ||
435 | fHeader->fNClus = fCaloClusters->GetEntriesFast(); | |
436 | fHeader->fTrackMult = trackMult; | |
b8b781cb | 437 | |
33d8da67 | 438 | fMCEvent = MCEvent(); |
80c98845 | 439 | if(fMCEvent){ |
33d8da67 | 440 | fStack = (AliStack*)fMCEvent->Stack(); |
4b6ab175 | 441 | if(this->fStack) |
442 | fHeader->fNMcParts = this->fStack->GetNtrack(); | |
7a7744db | 443 | else{ |
444 | printf("Stack not found\n"); | |
445 | fHeader->fNMcParts = 0; | |
446 | fAODMCParticles = (TClonesArray*)fVev->FindListObject("mcparticles"); | |
447 | } | |
448 | if(fAODMCParticles) | |
449 | fHeader->fNMcParts = fAODMCParticles->GetEntriesFast(); | |
80c98845 | 450 | } |
c5f3236f | 451 | else{ |
452 | if(fIsMC){ | |
453 | printf("ERROR: MC Event not available, returning...\n"); | |
454 | return; | |
455 | } | |
456 | } | |
33d8da67 | 457 | |
458 | ||
459 | FindConversions(); | |
4b6ab175 | 460 | if(this->fDebug) |
7a7744db | 461 | printf("FindConversions done\n"); |
33d8da67 | 462 | FillMyCells(); |
4b6ab175 | 463 | if(this->fDebug) |
7a7744db | 464 | printf("FillMyCells done\n"); |
33d8da67 | 465 | FillMyClusters(); |
4b6ab175 | 466 | if(this->fDebug) |
7a7744db | 467 | printf("FillMyClusters done\n"); |
33d8da67 | 468 | if(fCaloClustersNew) |
469 | FillMyAltClusters(); | |
3b37c011 | 470 | FillIsoTracks(); |
33d8da67 | 471 | if(fIsMC) |
472 | GetMcParts(); | |
c5f3236f | 473 | |
d1faab8e | 474 | if(this->fDebug && fIsMC) |
c5f3236f | 475 | printf("fMyMcParts nentries=%d",fMyMcParts->GetEntries()); |
33d8da67 | 476 | |
477 | fTree->Fill(); | |
478 | fSelTracks->Clear(); | |
479 | fSelPrimTracks->Clear(); | |
480 | fPhotConvArray->Clear(); | |
481 | fMyClusts->Clear(); | |
482 | if(fMyAltClusts) | |
483 | fMyAltClusts->Clear(); | |
484 | fMyCells->Clear(); | |
485 | fMyTracks->Clear(); | |
486 | if(fMyMcParts) | |
487 | fMyMcParts->Clear(); | |
c5f3236f | 488 | fMyMcIndex = 0; |
33d8da67 | 489 | fCaloClusters->Clear(); |
490 | if(fCaloClustersNew) | |
491 | fCaloClustersNew->Clear(); | |
7a7744db | 492 | if(fVCells) |
493 | fVCells = 0; | |
494 | // Post output data. | |
33d8da67 | 495 | PostData(1, fOutputList); |
735294e8 | 496 | PostData(2, fTree); |
33d8da67 | 497 | } |
498 | ||
499 | //________________________________________________________________________ | |
7a7744db | 500 | void AliAnalysisTaskEMCALPhoton::FindConversions() //WARNING: not ready to use with AODs |
33d8da67 | 501 | { |
f1c66719 | 502 | // Find conversion. |
7a7744db | 503 | if(!fESD)//not working with AODs yet |
504 | return; | |
4b6ab175 | 505 | if(this->fDebug) |
c5f3236f | 506 | printf("::FindConversions() starting\n"); |
33d8da67 | 507 | if(!fTrCuts) |
508 | return; | |
509 | Int_t iconv = 0; | |
7a7744db | 510 | Int_t nV0Orig = 0; |
511 | if(fESD) | |
512 | nV0Orig = fESD->GetNumberOfV0s(); | |
513 | if(fAOD) | |
514 | nV0Orig = fAOD->GetNumberOfV0s(); | |
4a4936e3 | 515 | Int_t nV0s = nV0Orig; |
7a7744db | 516 | Int_t ntracks = fSelTracks->GetEntriesFast(); |
517 | if(fRedoV0 && !fAOD){ | |
4a4936e3 | 518 | fESD->ResetV0s(); |
519 | AliV0vertexer lV0vtxer; | |
520 | lV0vtxer.Tracks2V0vertices(fESD); | |
521 | nV0s = fESD->GetNumberOfV0s(); | |
522 | } | |
33d8da67 | 523 | fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s); |
524 | for(Int_t iv0=0; iv0<nV0s; iv0++){ | |
525 | AliESDv0 * v0 = fESD->GetV0(iv0); | |
526 | if(!v0) | |
527 | continue; | |
528 | Double_t vpos[3]; | |
529 | fInvMassV0->Fill(v0->GetEffMass()); | |
530 | v0->GetXYZ(vpos[0], vpos[1], vpos[2]); | |
531 | Int_t ipos = v0->GetPindex(); | |
532 | Int_t ineg = v0->GetNindex(); | |
533 | if(ipos<0 || ipos> ntracks) | |
534 | continue; | |
535 | if(ineg<0 || ineg> ntracks) | |
536 | continue; | |
7a7744db | 537 | AliESDtrack *pos = static_cast<AliESDtrack*>(fSelTracks->At(ipos)); |
538 | AliESDtrack *neg = static_cast<AliESDtrack*>(fSelTracks->At(ineg)); | |
33d8da67 | 539 | if(!pos || !neg) |
540 | continue; | |
33d8da67 | 541 | /*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65) |
542 | continue;*/ | |
543 | const AliExternalTrackParam * paramPos = v0->GetParamP() ; | |
544 | const AliExternalTrackParam * paramNeg = v0->GetParamN() ; | |
545 | if(!paramPos || !paramNeg) | |
546 | continue; | |
547 | if(pos->GetSign() <0){//change tracks | |
548 | pos=neg ; | |
549 | neg=fESD->GetTrack(v0->GetPindex()) ; | |
550 | paramPos=paramNeg ; | |
551 | paramNeg=v0->GetParamP() ; | |
552 | } | |
553 | AliKFParticle negKF(*paramNeg,11); | |
554 | AliKFParticle posKF(*paramPos,-11); | |
555 | AliKFParticle photKF(negKF,posKF) ; | |
33d8da67 | 556 | |
557 | if( neg->GetKinkIndex(0) > 0 || pos->GetKinkIndex(0) > 0) | |
558 | continue ; | |
33d8da67 | 559 | |
560 | if(pos->GetSign() == neg->GetSign()){ | |
561 | fInvMassV0SS->Fill(photKF.GetMass()); | |
562 | continue; | |
563 | } | |
564 | fInvMassV0KF->Fill(photKF.GetMass()); | |
565 | TLorentzVector photLV(photKF.GetPx(), photKF.GetPy(),photKF.GetPz(), photKF.GetE()); | |
566 | if(photLV.M()>0.05 || photLV.M()<0) | |
567 | continue; | |
568 | fConversionVtxXY->Fill(vpos[0], vpos[1]);//negKF.GetX(), negKF.GetY()); | |
569 | Double_t convPhi = TMath::ATan2(photLV.Py(),photLV.Px()); | |
570 | if(convPhi<0) | |
571 | convPhi+=TMath::TwoPi(); | |
572 | TVector3 vecpos(vpos); | |
573 | Double_t v0Phi = 0; | |
574 | if(vecpos.Perp()>0) | |
575 | vecpos.Phi(); | |
576 | if(v0Phi<0) | |
577 | v0Phi+=TMath::TwoPi(); | |
b2d49404 | 578 | AliPhotonConvObj *myconv = static_cast<AliPhotonConvObj*>(fPhotConvArray->New(iconv++)); |
33d8da67 | 579 | myconv->fPt = photLV.Pt(); |
580 | myconv->fEta = photLV.Eta(); | |
581 | myconv->fPhi = convPhi; | |
582 | myconv->fVR = vecpos.Perp(); | |
583 | if(vecpos.Perp()>0) | |
584 | myconv->fVEta = vecpos.Eta(); | |
585 | myconv->fVPhi = v0Phi; | |
586 | myconv->fMass = photLV.M(); | |
587 | myconv->fMcLabel = -3; //WARNING: include the correct labeling | |
588 | //negative daughter | |
589 | myconv->fNegPt = negKF.GetPt(); | |
590 | myconv->fNegEta = negKF.GetEta(); | |
591 | Double_t trackPhi = negKF.GetPhi(); | |
592 | if(trackPhi<0) | |
593 | trackPhi+=TMath::TwoPi(); | |
594 | myconv->fNegPhi = trackPhi; | |
595 | myconv->fNegDedx = neg->GetTPCsignal(); | |
596 | myconv->fNegMcLabel = neg->GetLabel(); | |
597 | //negative daughter | |
598 | myconv->fPosPt = posKF.GetPt(); | |
599 | myconv->fPosEta = posKF.GetEta(); | |
600 | trackPhi = posKF.GetPhi(); | |
601 | if(trackPhi<0) | |
602 | trackPhi+=TMath::TwoPi(); | |
603 | myconv->fPosPhi = trackPhi; | |
604 | myconv->fPosDedx = pos->GetTPCsignal(); | |
605 | myconv->fPosMcLabel = pos->GetLabel(); | |
606 | } | |
4b6ab175 | 607 | if(this->fDebug) |
c5f3236f | 608 | printf("::FindConversions() returning...\n\n"); |
33d8da67 | 609 | return; |
610 | } | |
f1c66719 | 611 | |
33d8da67 | 612 | //________________________________________________________________________ |
613 | void AliAnalysisTaskEMCALPhoton::FillMyCells() | |
614 | { | |
f1c66719 | 615 | // Fill cells. |
4b6ab175 | 616 | if(this->fDebug) |
c5f3236f | 617 | printf("::FillMyCells() starting\n"); |
f1c66719 | 618 | |
7a7744db | 619 | if (!fVCells) |
33d8da67 | 620 | return; |
7a7744db | 621 | Int_t ncells = fVCells->GetNumberOfCells(); |
30bd295d | 622 | Int_t mcel = 0, maxcelid=-1; |
623 | Double_t maxcellE = 0, maxcellEta=0, maxcellPhi=0; | |
33d8da67 | 624 | for(Int_t icell = 0; icell<ncells; icell++){ |
7a7744db | 625 | Int_t absID = TMath::Abs(fVCells->GetCellNumber(icell)); |
b2d49404 | 626 | AliPhotonCellObj *mycell = static_cast<AliPhotonCellObj*>(fMyCells->New(mcel++)); |
33d8da67 | 627 | Float_t eta=-1, phi=-1; |
628 | if(!fGeom){ | |
40b40be6 | 629 | std::cout<<"+++fGeom not found! MyCells branch will not be filled for this event!+++\n"; |
33d8da67 | 630 | return; |
631 | } | |
632 | if(!fGeom) | |
633 | return; | |
40b40be6 | 634 | /*if(!fIsMC)*/fGeom->EtaPhiFromIndex(absID,eta,phi); |
30bd295d | 635 | if(maxcellE<fVCells->GetCellAmplitude(absID)){ |
636 | maxcellE = fVCells->GetCellAmplitude(absID); | |
637 | maxcellEta = eta; | |
638 | maxcellPhi = phi; | |
639 | maxcelid = absID; | |
640 | } | |
33d8da67 | 641 | Float_t theta = 2*TMath::ATan(TMath::Exp(-eta)); |
642 | mycell->fAbsID = absID; | |
7a7744db | 643 | mycell->fE = fVCells->GetCellAmplitude(absID); |
644 | mycell->fEt = fVCells->GetCellAmplitude(absID)*TMath::Sin(theta); | |
33d8da67 | 645 | mycell->fEta = eta; |
646 | mycell->fPhi = phi; | |
7a7744db | 647 | mycell->fTime = fVCells->GetCellTime(absID); |
33d8da67 | 648 | } |
4b6ab175 | 649 | if(this->fDebug) |
c5f3236f | 650 | printf("::FillMyCells() returning...\n\n"); |
33d8da67 | 651 | } |
f1c66719 | 652 | |
33d8da67 | 653 | //________________________________________________________________________ |
654 | void AliAnalysisTaskEMCALPhoton::FillMyClusters() | |
655 | { | |
f1c66719 | 656 | // Fill clusters. |
4b6ab175 | 657 | if(this->fDebug) |
c5f3236f | 658 | printf("::FillMyClusters() starting\n"); |
f1c66719 | 659 | |
7a7744db | 660 | if (!fCaloClusters){ |
661 | printf("CaloClusters is empty!\n"); | |
33d8da67 | 662 | return; |
7a7744db | 663 | } |
33d8da67 | 664 | Int_t nclus = fCaloClusters->GetEntries(); |
7a7744db | 665 | if(0==nclus) |
666 | printf("CaloClusters has ZERO entries\n"); | |
30bd295d | 667 | Int_t mcl = 0, maxcelid=-1; |
668 | Double_t maxcellE=0, maxcellEtac=0,maxcellPhic=0; | |
33d8da67 | 669 | for(Int_t ic=0; ic < nclus; ic++){ |
7a7744db | 670 | AliVCluster *clus = static_cast<AliVCluster*>(fCaloClusters->At(ic)); |
33d8da67 | 671 | if(!clus) |
672 | continue; | |
673 | if(!clus->IsEMCAL()) | |
674 | continue; | |
675 | if(clus->E() < fClusThresh) | |
676 | continue; | |
58acb8ab | 677 | if(fDebug) |
678 | printf("cluster %d survived\n", ic); | |
33d8da67 | 679 | Float_t pos[3]; |
680 | clus->GetPosition(pos); | |
681 | TVector3 cpos(pos); | |
682 | TString cellsAbsID; | |
b2d49404 | 683 | AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++)); |
33d8da67 | 684 | myclus->fE = clus->E(); |
685 | myclus->fEt = clus->E()*TMath::Sin(cpos.Theta()); | |
686 | myclus->fR = cpos.Perp(); | |
687 | myclus->fEta = cpos.Eta(); | |
688 | myclus->fPhi = cpos.Phi(); | |
689 | if(cpos.Phi()<0){ | |
690 | myclus->fPhi+=TMath::TwoPi(); | |
691 | } | |
692 | myclus->fN = clus->GetNCells(); | |
693 | Short_t id = -1; | |
694 | myclus->fEmax = GetMaxCellEnergy( clus, id); | |
695 | myclus->fIdmax = id; | |
30bd295d | 696 | if(maxcellE < myclus->fEmax){ |
697 | maxcellE = myclus->fEmax; | |
698 | maxcelid = id; | |
699 | maxcellEtac = cpos.Eta(); | |
700 | maxcellPhic = cpos.Phi(); | |
701 | } | |
7a7744db | 702 | myclus->fTmax = fVCells->GetCellTime(id); |
33d8da67 | 703 | myclus->fEcross = GetCrossEnergy( clus, id); |
704 | myclus->fDisp = clus->GetDispersion(); | |
705 | myclus->fM20 = clus->GetM20(); | |
706 | myclus->fM02 = clus->GetM02(); | |
707 | myclus->fTrDEta = clus->GetTrackDz(); | |
708 | myclus->fTrDPhi = clus->GetTrackDx(); | |
709 | myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.); | |
710 | myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.); | |
711 | myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.); | |
712 | myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.); | |
713 | myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.); | |
714 | myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.); | |
715 | myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.); | |
716 | myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.); | |
717 | for(Int_t icell=0;icell<myclus->fN;icell++){ | |
718 | Int_t absID = clus->GetCellAbsId(icell); | |
719 | cellsAbsID.Append(Form("%d",absID)); | |
720 | cellsAbsID.Append(";"); | |
721 | } | |
722 | myclus->fCellsAbsId = cellsAbsID; | |
723 | myclus->fMcLabel = clus->GetLabel(); | |
724 | Int_t matchIndex = clus->GetTrackMatchedIndex(); | |
7a7744db | 725 | if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){ |
33d8da67 | 726 | myclus->fTrEp = -1; |
727 | continue; | |
728 | } | |
7a7744db | 729 | AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex)); |
38faf770 | 730 | if(!track) |
33d8da67 | 731 | continue; |
7a7744db | 732 | if(fESD){ |
733 | if(!fPrTrCuts) | |
734 | continue; | |
735 | if(!fPrTrCuts->IsSelected(track)) | |
736 | continue; | |
737 | } | |
38faf770 | 738 | |
33d8da67 | 739 | myclus->fTrEp = clus->E()/track->P(); |
38faf770 | 740 | myclus->fTrDedx = track->GetTPCsignal(); |
33d8da67 | 741 | } |
30bd295d | 742 | if(this->fDebug){ |
743 | printf(" ---===+++ Max Cell among clusters: id=%d, E=%1.2f, eta-clus=%1.2f, phi-clus=%1.2f\n",maxcelid,maxcellE,maxcellEtac,maxcellPhic); | |
c5f3236f | 744 | printf("::FillMyClusters() returning...\n\n"); |
30bd295d | 745 | } |
33d8da67 | 746 | |
747 | } | |
748 | //________________________________________________________________________ | |
749 | void AliAnalysisTaskEMCALPhoton::FillMyAltClusters() | |
750 | { | |
f1c66719 | 751 | // Fill clusters. |
4b6ab175 | 752 | if(this->fDebug) |
c5f3236f | 753 | printf("::FillMyAltClusters() starting\n"); |
f1c66719 | 754 | |
33d8da67 | 755 | if(!fCaloClustersNew) |
756 | return; | |
757 | Int_t nclus = fCaloClustersNew->GetEntries(); | |
758 | Int_t mcl = 0; | |
759 | for(Int_t ic=0; ic < nclus; ic++){ | |
7a7744db | 760 | AliVCluster *clus = static_cast<AliVCluster*>(fCaloClustersNew->At(ic)); |
33d8da67 | 761 | if(!clus) |
762 | continue; | |
763 | if(!clus->IsEMCAL()) | |
764 | continue; | |
765 | if(clus->E() < fClusThresh) | |
766 | continue; | |
767 | Float_t pos[3]; | |
768 | clus->GetPosition(pos); | |
769 | TVector3 cpos(pos); | |
770 | TString cellsAbsID; | |
b2d49404 | 771 | AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++)); |
33d8da67 | 772 | myclus->fE = clus->E(); |
773 | myclus->fEt = clus->E()*TMath::Sin(cpos.Theta()); | |
774 | myclus->fR = cpos.Perp(); | |
775 | myclus->fEta = cpos.Eta(); | |
776 | myclus->fPhi = cpos.Phi(); | |
777 | if(cpos.Phi()<0){ | |
778 | myclus->fPhi+=TMath::TwoPi(); | |
779 | } | |
780 | myclus->fN = clus->GetNCells(); | |
781 | myclus->fDisp = clus->GetDispersion(); | |
782 | myclus->fM20 = clus->GetM20(); | |
783 | myclus->fM02 = clus->GetM02(); | |
784 | myclus->fTrDEta = clus->GetTrackDz(); | |
785 | myclus->fTrDPhi = clus->GetTrackDx(); | |
786 | myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.); | |
787 | myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.); | |
788 | myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.); | |
789 | myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.); | |
790 | for(Int_t icell=0;icell<myclus->fN;icell++){ | |
791 | Int_t absID = clus->GetCellAbsId(icell); | |
792 | cellsAbsID.Append(Form("%d",absID)); | |
793 | cellsAbsID.Append(";"); | |
794 | } | |
795 | myclus->fCellsAbsId = cellsAbsID; | |
796 | myclus->fMcLabel = clus->GetLabel(); | |
797 | Int_t matchIndex = clus->GetTrackMatchedIndex(); | |
7a7744db | 798 | if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){ |
33d8da67 | 799 | myclus->fTrEp = -1; |
800 | continue; | |
801 | } | |
7a7744db | 802 | AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex)); |
33d8da67 | 803 | if(!track){ |
804 | myclus->fTrEp = -1; | |
805 | continue; | |
806 | } | |
807 | if(!fPrTrCuts){ | |
808 | myclus->fTrEp = -1; | |
809 | continue; | |
810 | } | |
811 | if(!fPrTrCuts->IsSelected(track)){ | |
812 | myclus->fTrEp = -1; | |
813 | continue; | |
814 | } | |
815 | myclus->fTrEp = clus->E()/track->P(); | |
816 | } | |
4b6ab175 | 817 | if(this->fDebug) |
c5f3236f | 818 | printf("::FillMyAltClusters() returning...\n\n"); |
33d8da67 | 819 | |
820 | } | |
821 | //________________________________________________________________________ | |
3b37c011 | 822 | void AliAnalysisTaskEMCALPhoton::FillIsoTracks() |
33d8da67 | 823 | { |
f1c66719 | 824 | // Fill high pt tracks. |
4b6ab175 | 825 | if(this->fDebug) |
c5f3236f | 826 | printf("::FillIsoTracks() starting\n"); |
f1c66719 | 827 | |
33d8da67 | 828 | if(!fSelPrimTracks) |
829 | return; | |
830 | Int_t ntracks = fSelPrimTracks->GetEntries(); | |
831 | Int_t imtr = 0; | |
832 | for(Int_t it=0;it<ntracks; it++){ | |
7a7744db | 833 | AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(it)); |
33d8da67 | 834 | if(!track) |
835 | continue; | |
01166399 | 836 | /*if(track->Phi()<1.0 || track->Phi()>3.55) |
837 | continue;*/ | |
3b37c011 | 838 | if(TMath::Abs(track->Eta())>1.1) |
839 | continue; | |
a72860a2 | 840 | /*if(track->Pt()<3) |
841 | continue;*/ | |
842 | AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++)); | |
33d8da67 | 843 | mtr->fPt = track->Pt(); |
844 | mtr->fEta = track->Eta(); | |
845 | mtr->fPhi = track->Phi(); | |
3b37c011 | 846 | mtr->fCharge = track->Charge(); |
33d8da67 | 847 | mtr->fDedx = track->GetTPCsignal(); |
848 | mtr->fMcLabel = track->GetLabel(); | |
849 | } | |
4b6ab175 | 850 | if(this->fDebug) |
c5f3236f | 851 | printf("::FillIsoTracks() returning...\n\n"); |
33d8da67 | 852 | } |
f1c66719 | 853 | |
33d8da67 | 854 | //________________________________________________________________________ |
855 | void AliAnalysisTaskEMCALPhoton::GetMcParts() | |
856 | { | |
c5f3236f | 857 | // Get MC particles. |
4b6ab175 | 858 | if(this->fDebug) |
c5f3236f | 859 | printf("::GetMcParts() starting\n"); |
f1c66719 | 860 | |
4b6ab175 | 861 | if (!this->fStack && !fAODMCParticles) |
33d8da67 | 862 | return; |
4b6ab175 | 863 | if(this->fDebug) |
7a7744db | 864 | printf("either stack or aodmcpaticles exists\n"); |
865 | const AliVVertex *evtVtx = 0; | |
4b6ab175 | 866 | if(this->fStack) |
7a7744db | 867 | evtVtx = fMCEvent->GetPrimaryVertex(); |
868 | else | |
869 | printf("no such thing as mc vvtx\n"); | |
870 | /*if (!evtVtx) | |
871 | return;*/ | |
4b6ab175 | 872 | if(this->fDebug) |
7a7744db | 873 | printf("mc vvertex ok\n"); |
874 | Int_t nTracks = 0; | |
4b6ab175 | 875 | if(this->fStack) |
876 | nTracks = this->fStack->GetNtrack(); | |
7a7744db | 877 | else if(fAODMCParticles) |
878 | nTracks = fAODMCParticles->GetEntriesFast(); | |
879 | TParticle *mcP = 0; | |
880 | AliAODMCParticle *amcP = 0; | |
4b6ab175 | 881 | if(this->fDebug) |
7a7744db | 882 | printf("loop in the %d mc particles starting\n",nTracks); |
33d8da67 | 883 | for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) { |
4b6ab175 | 884 | if(this->fStack) |
885 | mcP = dynamic_cast<TParticle*>(this->fStack->Particle(iTrack)); | |
7a7744db | 886 | if(fAODMCParticles) |
887 | amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack)); | |
c5f3236f | 888 | |
33d8da67 | 889 | // primary particle |
7a7744db | 890 | Double_t dR = 0; |
891 | if(mcP){ | |
892 | dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) + | |
893 | (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) + | |
894 | (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ())); | |
895 | } | |
896 | if((dR > 0.5)) | |
33d8da67 | 897 | continue; |
c5f3236f | 898 | |
33d8da67 | 899 | |
900 | // kinematic cuts | |
7a7744db | 901 | Double_t pt = 0; |
902 | Double_t eta = 0; | |
903 | Double_t phi = 0; | |
904 | Int_t mother = -1; | |
905 | Int_t pdgcode = 0; | |
906 | if(mcP){ | |
907 | pt = mcP->Pt() ; | |
908 | eta = mcP->Eta(); | |
909 | phi = mcP->Phi(); | |
910 | mother = mcP->GetMother(0); | |
911 | pdgcode = mcP->GetPdgCode(); | |
912 | } else { | |
913 | if(amcP){ | |
914 | pt = amcP->Pt(); | |
915 | eta = amcP->Eta(); | |
916 | phi = amcP->Phi(); | |
917 | mother = amcP->GetMother(); | |
918 | pdgcode = amcP->GetPdgCode(); | |
919 | } else | |
920 | continue; | |
921 | } | |
c5f3236f | 922 | if (pt<0.5) |
33d8da67 | 923 | continue; |
c5f3236f | 924 | |
c5f3236f | 925 | if (TMath::Abs(eta)>0.7) |
33d8da67 | 926 | continue; |
c5f3236f | 927 | |
c5f3236f | 928 | if (phi<1.0||phi>3.3) |
33d8da67 | 929 | continue; |
c5f3236f | 930 | |
7a7744db | 931 | if (mother!=6 && mother!=7 ) |
c5f3236f | 932 | continue; |
933 | ||
934 | ||
33d8da67 | 935 | // pion or eta meson or direct photon |
7a7744db | 936 | if(pdgcode == 111) { |
937 | } else if(pdgcode == 221) { | |
938 | } else if(pdgcode == 22 ) { | |
80c98845 | 939 | } else { |
33d8da67 | 940 | continue; |
80c98845 | 941 | } |
c5f3236f | 942 | |
7a7744db | 943 | FillMcPart( fMyMcIndex++, iTrack); |
33d8da67 | 944 | } |
4b6ab175 | 945 | if(this->fDebug) |
c5f3236f | 946 | printf("::GetMcParts() returning...\n\n"); |
33d8da67 | 947 | } |
f1c66719 | 948 | |
33d8da67 | 949 | //________________________________________________________________________ |
7a7744db | 950 | void AliAnalysisTaskEMCALPhoton::FillMcPart( Int_t itrack, Int_t label) |
33d8da67 | 951 | { |
f1c66719 | 952 | // Fill MC particles. |
7a7744db | 953 | Int_t nTracks = 0; |
954 | TParticle *mcP = 0; | |
955 | AliAODMCParticle *amcP= 0; | |
4b6ab175 | 956 | if(this->fStack){ |
957 | nTracks = this->fStack->GetNtrack(); | |
958 | mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack)); | |
7a7744db | 959 | } |
960 | else if(fAODMCParticles){ | |
961 | nTracks = fAODMCParticles->GetEntriesFast(); | |
962 | amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(itrack)); | |
963 | } | |
4b6ab175 | 964 | if(this->fDebug) |
c5f3236f | 965 | printf("\t::FillMcParts() starting with label %d\n", itrack); |
7a7744db | 966 | TVector3 vmcv; |
967 | if(mcP) | |
968 | vmcv.SetXYZ(mcP->Vx(),mcP->Vy(), mcP->Vz()); | |
969 | else { | |
970 | if(amcP) | |
971 | vmcv.SetXYZ(amcP->Xv(),amcP->Yv(), amcP->Zv()); | |
972 | else | |
973 | return; | |
974 | } | |
975 | ||
80c98845 | 976 | AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(itrack)); |
c5f3236f | 977 | mcp->fLabel = label ; |
7a7744db | 978 | if(mcP){ |
979 | mcp->fPdg = mcP->GetPdgCode() ; | |
980 | mcp->fPt = mcP->Pt() ; | |
981 | mcp->fEta = mcP->Eta() ; | |
982 | mcp->fPhi = mcP->Phi() ; | |
983 | mcp->fMother = mcP->GetMother(0) ; | |
984 | mcp->fFirstD = mcP->GetFirstDaughter() ; | |
985 | mcp->fLastD = mcP->GetLastDaughter() ; | |
986 | mcp->fStatus = mcP->GetStatusCode(); | |
987 | } | |
988 | if(amcP){ | |
989 | mcp->fPdg = amcP->GetPdgCode() ; | |
990 | mcp->fPt = amcP->Pt() ; | |
991 | mcp->fEta = amcP->Eta() ; | |
992 | mcp->fPhi = amcP->Phi() ; | |
993 | mcp->fMother = amcP->GetMother() ; | |
994 | mcp->fFirstD = amcP->GetDaughter(0) ; | |
995 | mcp->fLastD = amcP->GetDaughter(amcP->GetNDaughters()-1) ; | |
996 | mcp->fStatus = amcP->GetStatus(); | |
997 | } | |
33d8da67 | 998 | mcp->fVR = vmcv.Perp(); |
999 | if(vmcv.Perp()>0){ | |
1000 | mcp->fVEta = vmcv.Eta() ; | |
1001 | mcp->fVPhi = vmcv.Phi() ; | |
1002 | } | |
7a7744db | 1003 | if(itrack == 8){ |
1004 | mcp->fIso = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.4 , 0.2); | |
1005 | mcp->fIso3 = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.3 , 0.2); | |
1006 | } | |
4b6ab175 | 1007 | if(this->fDebug) |
7a7744db | 1008 | printf("\t::FillMcParts(): label=%d, pdg=%d, pt=%1.1f, mom=%d, 1stD=%d,last=%d\n\t::FillMcParts() returning...\n\n", mcp->fLabel,mcp->fPdg,mcp->fPt,mcp->fMother,mcp->fFirstD,mcp->fLastD); |
1009 | for(Int_t id=mcp->fFirstD; id < mcp->fLastD; id++){ | |
1010 | if(id<=mcp->fMother) | |
c5f3236f | 1011 | continue; |
1012 | if(id<0 || id>nTracks) | |
1013 | continue; | |
7a7744db | 1014 | FillMcPart( fMyMcIndex++, id); |
c5f3236f | 1015 | } |
1016 | ||
80c98845 | 1017 | } |
1018 | //________________________________________________________________________ | |
7a7744db | 1019 | Double_t AliAnalysisTaskEMCALPhoton::GetMcIsolation( Int_t itrack, Double_t radius, Double_t ptcut) const |
80c98845 | 1020 | { |
4b6ab175 | 1021 | if(this->fDebug){ |
c5f3236f | 1022 | printf("\t\t::GetMcIsolation() starting\n"); |
7a7744db | 1023 | //printf("\t\t incoming particle: PDG = %d, itrack=%d;\n",mcP->GetPdgCode(),itrack); |
c5f3236f | 1024 | } |
6240502f | 1025 | if (!this->fStack && !this->fAODMCParticles && this->fDebug){ |
7a7744db | 1026 | printf("\t\t\tNo MC stack/array!\n"); |
80c98845 | 1027 | return -1; |
7a7744db | 1028 | } |
80c98845 | 1029 | if(itrack<6 || itrack>8) |
1030 | return -1; | |
4b6ab175 | 1031 | if(this->fDebug) |
7a7744db | 1032 | printf("\t\t\tparticle of interest selected\n"); |
1033 | TParticle *mcP = 0; | |
1034 | AliAODMCParticle *amcP = 0; | |
1035 | Int_t pdgcode = 0; | |
1036 | Float_t eta = 0; | |
1037 | Float_t phi = 0; | |
80c98845 | 1038 | Double_t sumpt=0; |
80c98845 | 1039 | Float_t dR; |
7a7744db | 1040 | Int_t nparts = 0; |
4b6ab175 | 1041 | if(this->fStack){ |
7a7744db | 1042 | nparts = fStack->GetNtrack(); |
4b6ab175 | 1043 | mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack)); |
7a7744db | 1044 | eta = mcP->Eta(); |
1045 | phi = mcP->Phi(); | |
1046 | pdgcode = mcP->GetPdgCode(); | |
1047 | } | |
4b6ab175 | 1048 | if(this->fAODMCParticles){ |
7a7744db | 1049 | nparts = fAODMCParticles->GetEntriesFast(); |
4b6ab175 | 1050 | amcP = dynamic_cast<AliAODMCParticle*>(this->fAODMCParticles->At(itrack)); |
1051 | if(amcP){ | |
1052 | eta = amcP->Eta(); | |
1053 | phi = amcP->Phi(); | |
1054 | pdgcode = amcP->GetPdgCode(); | |
1055 | } | |
7a7744db | 1056 | } |
1057 | if(pdgcode!=22) | |
1058 | return -1; | |
1059 | TParticle *mcisop = 0; | |
1060 | AliAODMCParticle *amcisop = 0; | |
80c98845 | 1061 | for(Int_t ip = 0; ip<nparts; ip++){ |
80c98845 | 1062 | if(ip==itrack) |
1063 | continue; | |
4b6ab175 | 1064 | if(this->fStack) |
1065 | mcisop = dynamic_cast<TParticle*>(this->fStack->Particle(ip)); | |
7a7744db | 1066 | if(fAODMCParticles) |
1067 | amcisop = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(ip)); | |
1068 | Int_t status = 0; | |
1069 | Int_t mother = 0; | |
1070 | Float_t pt = 0; | |
1071 | Float_t isophi = -99; | |
1072 | Float_t isoeta = -99; | |
1073 | TVector3 vmcv; | |
1074 | if(mcisop){ | |
1075 | status = mcisop->GetStatusCode(); | |
1076 | mother = mcisop->GetMother(0); | |
1077 | pt = mcisop->Pt(); | |
1078 | isophi = mcisop->Phi(); | |
1079 | isoeta = mcisop->Eta(); | |
1080 | vmcv.SetXYZ(mcisop->Vx(),mcisop->Vy(), mcisop->Vz()); | |
1081 | } | |
1082 | else { | |
1083 | if(amcisop){ | |
1084 | status = amcisop->GetStatus(); | |
1085 | mother = amcisop->GetMother(); | |
1086 | pt = amcisop->Pt(); | |
1087 | isophi = amcisop->Phi(); | |
1088 | isoeta = amcisop->Eta(); | |
1089 | vmcv.SetXYZ(amcisop->Xv(),amcisop->Yv(), amcisop->Zv()); | |
1090 | } | |
1091 | else | |
1092 | continue; | |
1093 | } | |
1094 | if(status!=1) | |
80c98845 | 1095 | continue; |
7a7744db | 1096 | if(mother == itrack) |
c5f3236f | 1097 | continue; |
7a7744db | 1098 | if(pt<ptcut) |
80c98845 | 1099 | continue; |
80c98845 | 1100 | if(vmcv.Perp()>1) |
1101 | continue; | |
7a7744db | 1102 | dR = TMath::Sqrt((phi-isophi)*(phi-isophi)+(eta-isoeta)*(eta-isoeta)); |
80c98845 | 1103 | if(dR>radius) |
1104 | continue; | |
7a7744db | 1105 | sumpt += pt; |
80c98845 | 1106 | } |
4b6ab175 | 1107 | if(this->fDebug) |
c5f3236f | 1108 | printf("\t\t::GetMcIsolation() returning value %f ...\n\n",sumpt); |
80c98845 | 1109 | return sumpt; |
c5f3236f | 1110 | } |
f1c66719 | 1111 | |
33d8da67 | 1112 | //________________________________________________________________________ |
1113 | Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const | |
1114 | { | |
1115 | // Compute isolation based on tracks. | |
4b6ab175 | 1116 | if(this->fDebug) |
c5f3236f | 1117 | printf("\t::GetTrackIsolation() starting\n"); |
1118 | ||
33d8da67 | 1119 | Double_t trkIsolation = 0; |
1120 | Double_t rad2 = radius*radius; | |
1121 | Int_t ntrks = fSelPrimTracks->GetEntries(); | |
1122 | for(Int_t j = 0; j<ntrks; ++j) { | |
1123 | AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j)); | |
1124 | if (!track) | |
1125 | continue; | |
1126 | if (track->Pt()<pt) | |
1127 | continue; | |
1128 | ||
1129 | Float_t eta = track->Eta(); | |
1130 | Float_t phi = track->Phi(); | |
1131 | Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi); | |
1132 | Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff; | |
1133 | if(dist>rad2) | |
1134 | continue; | |
1135 | trkIsolation += track->Pt(); | |
1136 | } | |
4b6ab175 | 1137 | if(this->fDebug) |
c5f3236f | 1138 | printf("\t::GetTrackIsolation() returning\n\n"); |
33d8da67 | 1139 | return trkIsolation; |
1140 | } | |
f1c66719 | 1141 | |
33d8da67 | 1142 | //________________________________________________________________________ |
1143 | Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const | |
1144 | { | |
f1c66719 | 1145 | // Get phi band. |
1146 | ||
33d8da67 | 1147 | if(!fSelPrimTracks) |
1148 | return 0; | |
1149 | Int_t ntracks = fSelPrimTracks->GetEntries(); | |
1150 | Double_t loweta = eta - R; | |
1151 | Double_t upeta = eta + R; | |
1152 | Double_t upphi = phi + R; | |
1153 | Double_t lowphi = phi - R; | |
1154 | Double_t et = 0; | |
1155 | for(Int_t itr=0; itr<ntracks; itr++){ | |
7a7744db | 1156 | AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(itr)); |
33d8da67 | 1157 | if(!track) |
1158 | continue; | |
1159 | if(track->Pt()<minpt) | |
1160 | continue; | |
1161 | if((track->Eta() < upeta) && (track->Eta() > loweta)) | |
1162 | continue; | |
1163 | if((track->Phi() > upphi) || (track->Phi() < lowphi)) | |
1164 | continue; | |
1165 | et+=track->Pt(); | |
1166 | } | |
1167 | return et; | |
1168 | } | |
f1c66719 | 1169 | |
33d8da67 | 1170 | //________________________________________________________________________ |
1171 | Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax) | |
1172 | { | |
1173 | // Calculate the energy of cross cells around the leading cell. | |
7a7744db | 1174 | if(!fVCells) |
33d8da67 | 1175 | return 0; |
1176 | ||
33d8da67 | 1177 | if (!fGeom) |
1178 | return 0; | |
1179 | ||
1180 | Int_t iSupMod = -1; | |
1181 | Int_t iTower = -1; | |
1182 | Int_t iIphi = -1; | |
1183 | Int_t iIeta = -1; | |
1184 | Int_t iphi = -1; | |
1185 | Int_t ieta = -1; | |
1186 | Int_t iphis = -1; | |
1187 | Int_t ietas = -1; | |
1188 | ||
1189 | Double_t crossEnergy = 0; | |
1190 | ||
1191 | fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta); | |
1192 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas); | |
1193 | ||
1194 | Int_t ncells = cluster->GetNCells(); | |
1195 | for (Int_t i=0; i<ncells; i++) { | |
1196 | Int_t cellAbsId = cluster->GetCellAbsId(i); | |
1197 | fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); | |
1198 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta); | |
1199 | Int_t aphidiff = TMath::Abs(iphi-iphis); | |
1200 | if (aphidiff>1) | |
1201 | continue; | |
1202 | Int_t aetadiff = TMath::Abs(ieta-ietas); | |
1203 | if (aetadiff>1) | |
1204 | continue; | |
1205 | if ( (aphidiff==1 && aetadiff==0) || | |
1206 | (aphidiff==0 && aetadiff==1) ) { | |
7a7744db | 1207 | crossEnergy += fVCells->GetCellAmplitude(cellAbsId); |
33d8da67 | 1208 | } |
1209 | } | |
1210 | ||
1211 | return crossEnergy; | |
1212 | } | |
1213 | ||
33d8da67 | 1214 | //________________________________________________________________________ |
1215 | Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const | |
1216 | { | |
1217 | // Get maximum energy of attached cell. | |
1218 | ||
1219 | id = -1; | |
7a7744db | 1220 | if(!fVCells) |
33d8da67 | 1221 | return 0; |
1222 | ||
1223 | Double_t maxe = 0; | |
1224 | Int_t ncells = cluster->GetNCells(); | |
1225 | for (Int_t i=0; i<ncells; i++) { | |
7a7744db | 1226 | Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); |
33d8da67 | 1227 | if (e>maxe) { |
1228 | maxe = e; | |
1229 | id = cluster->GetCellAbsId(i); | |
1230 | } | |
1231 | } | |
1232 | return maxe; | |
1233 | } | |
f1c66719 | 1234 | |
33d8da67 | 1235 | //________________________________________________________________________ |
1236 | void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *) | |
1237 | { | |
33d8da67 | 1238 | // Called once at the end of the query |
ec681cf3 | 1239 | /* if(fIsGrid) |
1240 | return;*/ | |
f1c66719 | 1241 | if (fTree) { |
c305511b | 1242 | printf("***tree %s being saved***\n",fTree->GetName()); |
55dcd0ae | 1243 | TFile *f = OpenFile(2); |
33d8da67 | 1244 | TDirectory::TContext context(f); |
1245 | if (f) | |
1246 | fTree->Write(); | |
1247 | } | |
33d8da67 | 1248 | } |