]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhoton.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Correction in TObjArray* AliTPythia8::ImportParticles(Option_t* /* option */)
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALPhoton.cxx
CommitLineData
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 52using std::cout;
53using std::endl;
54
33d8da67 55ClassImp(AliAnalysisTaskEMCALPhoton)
56
57//________________________________________________________________________
f1c66719 58AliAnalysisTaskEMCALPhoton::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//________________________________________________________________________
110AliAnalysisTaskEMCALPhoton::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//________________________________________________________________________
169void 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//________________________________________________________________________
284void 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();
7a7744db 387 fHeader->fTrClassMask = fVev->GetHeader()->GetTriggerMask();
388 fHeader->fTrCluster = fVev->GetHeader()->GetTriggerCluster();
33d8da67 389 AliCentrality *cent = InputEvent()->GetCentrality();
390 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
391 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
392 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
7a7744db 393
5b19b2cc 394 if(fDebug)
395 printf("AliEMCALgeo file found\n");
396 TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices");
7a7744db 397 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
398 if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
399 break;
400 /*if(fESD)
33d8da67 401 fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
7a7744db 402 else*/
403 // if(event->GetEMCALMatrix(mod))
404 fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod);
405 fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod);
406 }
407 Int_t trackMult = 0;
408 if(fESD){
409 AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
410 trackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
4b6ab175 411 if(this->fDebug)
7a7744db 412 printf("ESD Track mult= %d\n",trackMult);
413 }
414 else if(fAOD){
415 trackMult = Ntracks;
4b6ab175 416 if(this->fDebug)
7a7744db 417 printf("AOD Track mult= %d\n",trackMult);
418 }
419 if(fESD){
420 TList *l = fESD->GetList();
421 fCaloClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
422 fVCells = fESD->GetEMCALCells();
423 fHeader->fNCells = fVCells->GetNumberOfCells();
4b6ab175 424 if(this->fDebug)
7a7744db 425 printf("ESD cluster mult= %d\n",fCaloClusters->GetEntriesFast());
426 }
427 else if(fAOD){
428 fCaloClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters());
429 fVCells = fAOD->GetEMCALCells();
430 fHeader->fNCells = fVCells->GetNumberOfCells();
4b6ab175 431 if(this->fDebug)
7a7744db 432 printf("AOD cluster mult= %d\n",fCaloClusters->GetEntriesFast());
33d8da67 433 }
7a7744db 434
435
436 fHeader->fNClus = fCaloClusters->GetEntriesFast();
437 fHeader->fTrackMult = trackMult;
b8b781cb 438
33d8da67 439 fMCEvent = MCEvent();
80c98845 440 if(fMCEvent){
33d8da67 441 fStack = (AliStack*)fMCEvent->Stack();
4b6ab175 442 if(this->fStack)
443 fHeader->fNMcParts = this->fStack->GetNtrack();
7a7744db 444 else{
445 printf("Stack not found\n");
446 fHeader->fNMcParts = 0;
447 fAODMCParticles = (TClonesArray*)fVev->FindListObject("mcparticles");
448 }
449 if(fAODMCParticles)
450 fHeader->fNMcParts = fAODMCParticles->GetEntriesFast();
80c98845 451 }
c5f3236f 452 else{
453 if(fIsMC){
454 printf("ERROR: MC Event not available, returning...\n");
455 return;
456 }
457 }
33d8da67 458
459
460 FindConversions();
4b6ab175 461 if(this->fDebug)
7a7744db 462 printf("FindConversions done\n");
33d8da67 463 FillMyCells();
4b6ab175 464 if(this->fDebug)
7a7744db 465 printf("FillMyCells done\n");
33d8da67 466 FillMyClusters();
4b6ab175 467 if(this->fDebug)
7a7744db 468 printf("FillMyClusters done\n");
33d8da67 469 if(fCaloClustersNew)
470 FillMyAltClusters();
3b37c011 471 FillIsoTracks();
33d8da67 472 if(fIsMC)
473 GetMcParts();
c5f3236f 474
d1faab8e 475 if(this->fDebug && fIsMC)
c5f3236f 476 printf("fMyMcParts nentries=%d",fMyMcParts->GetEntries());
33d8da67 477
478 fTree->Fill();
479 fSelTracks->Clear();
480 fSelPrimTracks->Clear();
481 fPhotConvArray->Clear();
482 fMyClusts->Clear();
483 if(fMyAltClusts)
484 fMyAltClusts->Clear();
485 fMyCells->Clear();
486 fMyTracks->Clear();
487 if(fMyMcParts)
488 fMyMcParts->Clear();
c5f3236f 489 fMyMcIndex = 0;
33d8da67 490 fCaloClusters->Clear();
491 if(fCaloClustersNew)
492 fCaloClustersNew->Clear();
7a7744db 493 if(fVCells)
494 fVCells = 0;
495 // Post output data.
33d8da67 496 PostData(1, fOutputList);
735294e8 497 PostData(2, fTree);
33d8da67 498}
499
500//________________________________________________________________________
7a7744db 501void AliAnalysisTaskEMCALPhoton::FindConversions() //WARNING: not ready to use with AODs
33d8da67 502{
f1c66719 503 // Find conversion.
7a7744db 504 if(!fESD)//not working with AODs yet
505 return;
4b6ab175 506 if(this->fDebug)
c5f3236f 507 printf("::FindConversions() starting\n");
33d8da67 508 if(!fTrCuts)
509 return;
510 Int_t iconv = 0;
7a7744db 511 Int_t nV0Orig = 0;
512 if(fESD)
513 nV0Orig = fESD->GetNumberOfV0s();
514 if(fAOD)
515 nV0Orig = fAOD->GetNumberOfV0s();
4a4936e3 516 Int_t nV0s = nV0Orig;
7a7744db 517 Int_t ntracks = fSelTracks->GetEntriesFast();
518 if(fRedoV0 && !fAOD){
4a4936e3 519 fESD->ResetV0s();
520 AliV0vertexer lV0vtxer;
521 lV0vtxer.Tracks2V0vertices(fESD);
522 nV0s = fESD->GetNumberOfV0s();
523 }
33d8da67 524 fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s);
525 for(Int_t iv0=0; iv0<nV0s; iv0++){
526 AliESDv0 * v0 = fESD->GetV0(iv0);
527 if(!v0)
528 continue;
529 Double_t vpos[3];
530 fInvMassV0->Fill(v0->GetEffMass());
531 v0->GetXYZ(vpos[0], vpos[1], vpos[2]);
532 Int_t ipos = v0->GetPindex();
533 Int_t ineg = v0->GetNindex();
534 if(ipos<0 || ipos> ntracks)
535 continue;
536 if(ineg<0 || ineg> ntracks)
537 continue;
7a7744db 538 AliESDtrack *pos = static_cast<AliESDtrack*>(fSelTracks->At(ipos));
539 AliESDtrack *neg = static_cast<AliESDtrack*>(fSelTracks->At(ineg));
33d8da67 540 if(!pos || !neg)
541 continue;
33d8da67 542 /*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65)
543 continue;*/
544 const AliExternalTrackParam * paramPos = v0->GetParamP() ;
545 const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
546 if(!paramPos || !paramNeg)
547 continue;
548 if(pos->GetSign() <0){//change tracks
549 pos=neg ;
550 neg=fESD->GetTrack(v0->GetPindex()) ;
551 paramPos=paramNeg ;
552 paramNeg=v0->GetParamP() ;
553 }
554 AliKFParticle negKF(*paramNeg,11);
555 AliKFParticle posKF(*paramPos,-11);
556 AliKFParticle photKF(negKF,posKF) ;
33d8da67 557
558 if( neg->GetKinkIndex(0) > 0 || pos->GetKinkIndex(0) > 0)
559 continue ;
33d8da67 560
561 if(pos->GetSign() == neg->GetSign()){
562 fInvMassV0SS->Fill(photKF.GetMass());
563 continue;
564 }
565 fInvMassV0KF->Fill(photKF.GetMass());
566 TLorentzVector photLV(photKF.GetPx(), photKF.GetPy(),photKF.GetPz(), photKF.GetE());
567 if(photLV.M()>0.05 || photLV.M()<0)
568 continue;
569 fConversionVtxXY->Fill(vpos[0], vpos[1]);//negKF.GetX(), negKF.GetY());
570 Double_t convPhi = TMath::ATan2(photLV.Py(),photLV.Px());
571 if(convPhi<0)
572 convPhi+=TMath::TwoPi();
573 TVector3 vecpos(vpos);
574 Double_t v0Phi = 0;
575 if(vecpos.Perp()>0)
576 vecpos.Phi();
577 if(v0Phi<0)
578 v0Phi+=TMath::TwoPi();
b2d49404 579 AliPhotonConvObj *myconv = static_cast<AliPhotonConvObj*>(fPhotConvArray->New(iconv++));
33d8da67 580 myconv->fPt = photLV.Pt();
581 myconv->fEta = photLV.Eta();
582 myconv->fPhi = convPhi;
583 myconv->fVR = vecpos.Perp();
584 if(vecpos.Perp()>0)
585 myconv->fVEta = vecpos.Eta();
586 myconv->fVPhi = v0Phi;
587 myconv->fMass = photLV.M();
588 myconv->fMcLabel = -3; //WARNING: include the correct labeling
589 //negative daughter
590 myconv->fNegPt = negKF.GetPt();
591 myconv->fNegEta = negKF.GetEta();
592 Double_t trackPhi = negKF.GetPhi();
593 if(trackPhi<0)
594 trackPhi+=TMath::TwoPi();
595 myconv->fNegPhi = trackPhi;
596 myconv->fNegDedx = neg->GetTPCsignal();
597 myconv->fNegMcLabel = neg->GetLabel();
598 //negative daughter
599 myconv->fPosPt = posKF.GetPt();
600 myconv->fPosEta = posKF.GetEta();
601 trackPhi = posKF.GetPhi();
602 if(trackPhi<0)
603 trackPhi+=TMath::TwoPi();
604 myconv->fPosPhi = trackPhi;
605 myconv->fPosDedx = pos->GetTPCsignal();
606 myconv->fPosMcLabel = pos->GetLabel();
607 }
4b6ab175 608 if(this->fDebug)
c5f3236f 609 printf("::FindConversions() returning...\n\n");
33d8da67 610 return;
611}
f1c66719 612
33d8da67 613//________________________________________________________________________
614void AliAnalysisTaskEMCALPhoton::FillMyCells()
615{
f1c66719 616 // Fill cells.
4b6ab175 617 if(this->fDebug)
c5f3236f 618 printf("::FillMyCells() starting\n");
f1c66719 619
7a7744db 620 if (!fVCells)
33d8da67 621 return;
7a7744db 622 Int_t ncells = fVCells->GetNumberOfCells();
33d8da67 623 Int_t mcel = 0;
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);
33d8da67 635 Float_t theta = 2*TMath::ATan(TMath::Exp(-eta));
636 mycell->fAbsID = absID;
7a7744db 637 mycell->fE = fVCells->GetCellAmplitude(absID);
638 mycell->fEt = fVCells->GetCellAmplitude(absID)*TMath::Sin(theta);
33d8da67 639 mycell->fEta = eta;
640 mycell->fPhi = phi;
7a7744db 641 mycell->fTime = fVCells->GetCellTime(absID);
33d8da67 642 }
4b6ab175 643 if(this->fDebug)
c5f3236f 644 printf("::FillMyCells() returning...\n\n");
33d8da67 645}
f1c66719 646
33d8da67 647//________________________________________________________________________
648void AliAnalysisTaskEMCALPhoton::FillMyClusters()
649{
f1c66719 650 // Fill clusters.
4b6ab175 651 if(this->fDebug)
c5f3236f 652 printf("::FillMyClusters() starting\n");
f1c66719 653
7a7744db 654 if (!fCaloClusters){
655 printf("CaloClusters is empty!\n");
33d8da67 656 return;
7a7744db 657 }
33d8da67 658 Int_t nclus = fCaloClusters->GetEntries();
7a7744db 659 if(0==nclus)
660 printf("CaloClusters has ZERO entries\n");
33d8da67 661 Int_t mcl = 0;
662 for(Int_t ic=0; ic < nclus; ic++){
7a7744db 663 AliVCluster *clus = static_cast<AliVCluster*>(fCaloClusters->At(ic));
33d8da67 664 if(!clus)
665 continue;
666 if(!clus->IsEMCAL())
667 continue;
668 if(clus->E() < fClusThresh)
669 continue;
670 Float_t pos[3];
671 clus->GetPosition(pos);
672 TVector3 cpos(pos);
673 TString cellsAbsID;
b2d49404 674 AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++));
33d8da67 675 myclus->fE = clus->E();
676 myclus->fEt = clus->E()*TMath::Sin(cpos.Theta());
677 myclus->fR = cpos.Perp();
678 myclus->fEta = cpos.Eta();
679 myclus->fPhi = cpos.Phi();
680 if(cpos.Phi()<0){
681 myclus->fPhi+=TMath::TwoPi();
682 }
683 myclus->fN = clus->GetNCells();
684 Short_t id = -1;
685 myclus->fEmax = GetMaxCellEnergy( clus, id);
686 myclus->fIdmax = id;
7a7744db 687 myclus->fTmax = fVCells->GetCellTime(id);
33d8da67 688 myclus->fEcross = GetCrossEnergy( clus, id);
689 myclus->fDisp = clus->GetDispersion();
690 myclus->fM20 = clus->GetM20();
691 myclus->fM02 = clus->GetM02();
692 myclus->fTrDEta = clus->GetTrackDz();
693 myclus->fTrDPhi = clus->GetTrackDx();
694 myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
695 myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
696 myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
697 myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
698 myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.);
699 myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.);
700 myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.);
701 myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.);
702 for(Int_t icell=0;icell<myclus->fN;icell++){
703 Int_t absID = clus->GetCellAbsId(icell);
704 cellsAbsID.Append(Form("%d",absID));
705 cellsAbsID.Append(";");
706 }
707 myclus->fCellsAbsId = cellsAbsID;
708 myclus->fMcLabel = clus->GetLabel();
709 Int_t matchIndex = clus->GetTrackMatchedIndex();
7a7744db 710 if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){
33d8da67 711 myclus->fTrEp = -1;
712 continue;
713 }
7a7744db 714 AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex));
38faf770 715 if(!track)
33d8da67 716 continue;
7a7744db 717 if(fESD){
718 if(!fPrTrCuts)
719 continue;
720 if(!fPrTrCuts->IsSelected(track))
721 continue;
722 }
38faf770 723
33d8da67 724 myclus->fTrEp = clus->E()/track->P();
38faf770 725 myclus->fTrDedx = track->GetTPCsignal();
33d8da67 726 }
4b6ab175 727 if(this->fDebug)
c5f3236f 728 printf("::FillMyClusters() returning...\n\n");
33d8da67 729
730}
731//________________________________________________________________________
732void AliAnalysisTaskEMCALPhoton::FillMyAltClusters()
733{
f1c66719 734 // Fill clusters.
4b6ab175 735 if(this->fDebug)
c5f3236f 736 printf("::FillMyAltClusters() starting\n");
f1c66719 737
33d8da67 738 if(!fCaloClustersNew)
739 return;
740 Int_t nclus = fCaloClustersNew->GetEntries();
741 Int_t mcl = 0;
742 for(Int_t ic=0; ic < nclus; ic++){
7a7744db 743 AliVCluster *clus = static_cast<AliVCluster*>(fCaloClustersNew->At(ic));
33d8da67 744 if(!clus)
745 continue;
746 if(!clus->IsEMCAL())
747 continue;
748 if(clus->E() < fClusThresh)
749 continue;
750 Float_t pos[3];
751 clus->GetPosition(pos);
752 TVector3 cpos(pos);
753 TString cellsAbsID;
b2d49404 754 AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++));
33d8da67 755 myclus->fE = clus->E();
756 myclus->fEt = clus->E()*TMath::Sin(cpos.Theta());
757 myclus->fR = cpos.Perp();
758 myclus->fEta = cpos.Eta();
759 myclus->fPhi = cpos.Phi();
760 if(cpos.Phi()<0){
761 myclus->fPhi+=TMath::TwoPi();
762 }
763 myclus->fN = clus->GetNCells();
764 myclus->fDisp = clus->GetDispersion();
765 myclus->fM20 = clus->GetM20();
766 myclus->fM02 = clus->GetM02();
767 myclus->fTrDEta = clus->GetTrackDz();
768 myclus->fTrDPhi = clus->GetTrackDx();
769 myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
770 myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
771 myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
772 myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
773 for(Int_t icell=0;icell<myclus->fN;icell++){
774 Int_t absID = clus->GetCellAbsId(icell);
775 cellsAbsID.Append(Form("%d",absID));
776 cellsAbsID.Append(";");
777 }
778 myclus->fCellsAbsId = cellsAbsID;
779 myclus->fMcLabel = clus->GetLabel();
780 Int_t matchIndex = clus->GetTrackMatchedIndex();
7a7744db 781 if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){
33d8da67 782 myclus->fTrEp = -1;
783 continue;
784 }
7a7744db 785 AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex));
33d8da67 786 if(!track){
787 myclus->fTrEp = -1;
788 continue;
789 }
790 if(!fPrTrCuts){
791 myclus->fTrEp = -1;
792 continue;
793 }
794 if(!fPrTrCuts->IsSelected(track)){
795 myclus->fTrEp = -1;
796 continue;
797 }
798 myclus->fTrEp = clus->E()/track->P();
799 }
4b6ab175 800 if(this->fDebug)
c5f3236f 801 printf("::FillMyAltClusters() returning...\n\n");
33d8da67 802
803}
804//________________________________________________________________________
3b37c011 805void AliAnalysisTaskEMCALPhoton::FillIsoTracks()
33d8da67 806{
f1c66719 807 // Fill high pt tracks.
4b6ab175 808 if(this->fDebug)
c5f3236f 809 printf("::FillIsoTracks() starting\n");
f1c66719 810
33d8da67 811 if(!fSelPrimTracks)
812 return;
813 Int_t ntracks = fSelPrimTracks->GetEntries();
814 Int_t imtr = 0;
815 for(Int_t it=0;it<ntracks; it++){
7a7744db 816 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(it));
33d8da67 817 if(!track)
818 continue;
01166399 819 /*if(track->Phi()<1.0 || track->Phi()>3.55)
820 continue;*/
3b37c011 821 if(TMath::Abs(track->Eta())>1.1)
822 continue;
a72860a2 823 /*if(track->Pt()<3)
824 continue;*/
825 AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
33d8da67 826 mtr->fPt = track->Pt();
827 mtr->fEta = track->Eta();
828 mtr->fPhi = track->Phi();
3b37c011 829 mtr->fCharge = track->Charge();
33d8da67 830 mtr->fDedx = track->GetTPCsignal();
831 mtr->fMcLabel = track->GetLabel();
832 }
4b6ab175 833 if(this->fDebug)
c5f3236f 834 printf("::FillIsoTracks() returning...\n\n");
33d8da67 835}
f1c66719 836
33d8da67 837//________________________________________________________________________
838void AliAnalysisTaskEMCALPhoton::GetMcParts()
839{
c5f3236f 840 // Get MC particles.
4b6ab175 841 if(this->fDebug)
c5f3236f 842 printf("::GetMcParts() starting\n");
f1c66719 843
4b6ab175 844 if (!this->fStack && !fAODMCParticles)
33d8da67 845 return;
4b6ab175 846 if(this->fDebug)
7a7744db 847 printf("either stack or aodmcpaticles exists\n");
848 const AliVVertex *evtVtx = 0;
4b6ab175 849 if(this->fStack)
7a7744db 850 evtVtx = fMCEvent->GetPrimaryVertex();
851 else
852 printf("no such thing as mc vvtx\n");
853 /*if (!evtVtx)
854 return;*/
4b6ab175 855 if(this->fDebug)
7a7744db 856 printf("mc vvertex ok\n");
857 Int_t nTracks = 0;
4b6ab175 858 if(this->fStack)
859 nTracks = this->fStack->GetNtrack();
7a7744db 860 else if(fAODMCParticles)
861 nTracks = fAODMCParticles->GetEntriesFast();
862 TParticle *mcP = 0;
863 AliAODMCParticle *amcP = 0;
4b6ab175 864 if(this->fDebug)
7a7744db 865 printf("loop in the %d mc particles starting\n",nTracks);
33d8da67 866 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
4b6ab175 867 if(this->fStack)
868 mcP = dynamic_cast<TParticle*>(this->fStack->Particle(iTrack));
7a7744db 869 if(fAODMCParticles)
870 amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
c5f3236f 871
33d8da67 872 // primary particle
7a7744db 873 Double_t dR = 0;
874 if(mcP){
875 dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
876 (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
877 (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
878 }
879 if((dR > 0.5))
33d8da67 880 continue;
c5f3236f 881
33d8da67 882
883 // kinematic cuts
7a7744db 884 Double_t pt = 0;
885 Double_t eta = 0;
886 Double_t phi = 0;
887 Int_t mother = -1;
888 Int_t pdgcode = 0;
889 if(mcP){
890 pt = mcP->Pt() ;
891 eta = mcP->Eta();
892 phi = mcP->Phi();
893 mother = mcP->GetMother(0);
894 pdgcode = mcP->GetPdgCode();
895 } else {
896 if(amcP){
897 pt = amcP->Pt();
898 eta = amcP->Eta();
899 phi = amcP->Phi();
900 mother = amcP->GetMother();
901 pdgcode = amcP->GetPdgCode();
902 } else
903 continue;
904 }
c5f3236f 905 if (pt<0.5)
33d8da67 906 continue;
c5f3236f 907
c5f3236f 908 if (TMath::Abs(eta)>0.7)
33d8da67 909 continue;
c5f3236f 910
c5f3236f 911 if (phi<1.0||phi>3.3)
33d8da67 912 continue;
c5f3236f 913
7a7744db 914 if (mother!=6 && mother!=7 )
c5f3236f 915 continue;
916
917
33d8da67 918 // pion or eta meson or direct photon
7a7744db 919 if(pdgcode == 111) {
920 } else if(pdgcode == 221) {
921 } else if(pdgcode == 22 ) {
80c98845 922 } else {
33d8da67 923 continue;
80c98845 924 }
c5f3236f 925
7a7744db 926 FillMcPart( fMyMcIndex++, iTrack);
33d8da67 927 }
4b6ab175 928 if(this->fDebug)
c5f3236f 929 printf("::GetMcParts() returning...\n\n");
33d8da67 930}
f1c66719 931
33d8da67 932//________________________________________________________________________
7a7744db 933void AliAnalysisTaskEMCALPhoton::FillMcPart( Int_t itrack, Int_t label)
33d8da67 934{
f1c66719 935 // Fill MC particles.
7a7744db 936 Int_t nTracks = 0;
937 TParticle *mcP = 0;
938 AliAODMCParticle *amcP= 0;
4b6ab175 939 if(this->fStack){
940 nTracks = this->fStack->GetNtrack();
941 mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack));
7a7744db 942 }
943 else if(fAODMCParticles){
944 nTracks = fAODMCParticles->GetEntriesFast();
945 amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(itrack));
946 }
4b6ab175 947 if(this->fDebug)
c5f3236f 948 printf("\t::FillMcParts() starting with label %d\n", itrack);
7a7744db 949 TVector3 vmcv;
950 if(mcP)
951 vmcv.SetXYZ(mcP->Vx(),mcP->Vy(), mcP->Vz());
952 else {
953 if(amcP)
954 vmcv.SetXYZ(amcP->Xv(),amcP->Yv(), amcP->Zv());
955 else
956 return;
957 }
958
80c98845 959 AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(itrack));
c5f3236f 960 mcp->fLabel = label ;
7a7744db 961 if(mcP){
962 mcp->fPdg = mcP->GetPdgCode() ;
963 mcp->fPt = mcP->Pt() ;
964 mcp->fEta = mcP->Eta() ;
965 mcp->fPhi = mcP->Phi() ;
966 mcp->fMother = mcP->GetMother(0) ;
967 mcp->fFirstD = mcP->GetFirstDaughter() ;
968 mcp->fLastD = mcP->GetLastDaughter() ;
969 mcp->fStatus = mcP->GetStatusCode();
970 }
971 if(amcP){
972 mcp->fPdg = amcP->GetPdgCode() ;
973 mcp->fPt = amcP->Pt() ;
974 mcp->fEta = amcP->Eta() ;
975 mcp->fPhi = amcP->Phi() ;
976 mcp->fMother = amcP->GetMother() ;
977 mcp->fFirstD = amcP->GetDaughter(0) ;
978 mcp->fLastD = amcP->GetDaughter(amcP->GetNDaughters()-1) ;
979 mcp->fStatus = amcP->GetStatus();
980 }
33d8da67 981 mcp->fVR = vmcv.Perp();
982 if(vmcv.Perp()>0){
983 mcp->fVEta = vmcv.Eta() ;
984 mcp->fVPhi = vmcv.Phi() ;
985 }
7a7744db 986 if(itrack == 8){
987 mcp->fIso = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.4 , 0.2);
988 mcp->fIso3 = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.3 , 0.2);
989 }
4b6ab175 990 if(this->fDebug)
7a7744db 991 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);
992 for(Int_t id=mcp->fFirstD; id < mcp->fLastD; id++){
993 if(id<=mcp->fMother)
c5f3236f 994 continue;
995 if(id<0 || id>nTracks)
996 continue;
7a7744db 997 FillMcPart( fMyMcIndex++, id);
c5f3236f 998 }
999
80c98845 1000}
1001//________________________________________________________________________
7a7744db 1002Double_t AliAnalysisTaskEMCALPhoton::GetMcIsolation( Int_t itrack, Double_t radius, Double_t ptcut) const
80c98845 1003{
4b6ab175 1004 if(this->fDebug){
c5f3236f 1005 printf("\t\t::GetMcIsolation() starting\n");
7a7744db 1006 //printf("\t\t incoming particle: PDG = %d, itrack=%d;\n",mcP->GetPdgCode(),itrack);
c5f3236f 1007 }
4b6ab175 1008 if (!this->fStack && !this->fAODMCParticles){
7a7744db 1009 printf("\t\t\tNo MC stack/array!\n");
80c98845 1010 return -1;
7a7744db 1011 }
80c98845 1012 if(itrack<6 || itrack>8)
1013 return -1;
4b6ab175 1014 if(this->fDebug)
7a7744db 1015 printf("\t\t\tparticle of interest selected\n");
1016 TParticle *mcP = 0;
1017 AliAODMCParticle *amcP = 0;
1018 Int_t pdgcode = 0;
1019 Float_t eta = 0;
1020 Float_t phi = 0;
80c98845 1021 Double_t sumpt=0;
80c98845 1022 Float_t dR;
7a7744db 1023 Int_t nparts = 0;
4b6ab175 1024 if(this->fStack){
7a7744db 1025 nparts = fStack->GetNtrack();
4b6ab175 1026 mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack));
7a7744db 1027 eta = mcP->Eta();
1028 phi = mcP->Phi();
1029 pdgcode = mcP->GetPdgCode();
1030 }
4b6ab175 1031 if(this->fAODMCParticles){
7a7744db 1032 nparts = fAODMCParticles->GetEntriesFast();
4b6ab175 1033 amcP = dynamic_cast<AliAODMCParticle*>(this->fAODMCParticles->At(itrack));
1034 if(amcP){
1035 eta = amcP->Eta();
1036 phi = amcP->Phi();
1037 pdgcode = amcP->GetPdgCode();
1038 }
7a7744db 1039 }
1040 if(pdgcode!=22)
1041 return -1;
1042 TParticle *mcisop = 0;
1043 AliAODMCParticle *amcisop = 0;
80c98845 1044 for(Int_t ip = 0; ip<nparts; ip++){
80c98845 1045 if(ip==itrack)
1046 continue;
4b6ab175 1047 if(this->fStack)
1048 mcisop = dynamic_cast<TParticle*>(this->fStack->Particle(ip));
7a7744db 1049 if(fAODMCParticles)
1050 amcisop = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(ip));
1051 Int_t status = 0;
1052 Int_t mother = 0;
1053 Float_t pt = 0;
1054 Float_t isophi = -99;
1055 Float_t isoeta = -99;
1056 TVector3 vmcv;
1057 if(mcisop){
1058 status = mcisop->GetStatusCode();
1059 mother = mcisop->GetMother(0);
1060 pt = mcisop->Pt();
1061 isophi = mcisop->Phi();
1062 isoeta = mcisop->Eta();
1063 vmcv.SetXYZ(mcisop->Vx(),mcisop->Vy(), mcisop->Vz());
1064 }
1065 else {
1066 if(amcisop){
1067 status = amcisop->GetStatus();
1068 mother = amcisop->GetMother();
1069 pt = amcisop->Pt();
1070 isophi = amcisop->Phi();
1071 isoeta = amcisop->Eta();
1072 vmcv.SetXYZ(amcisop->Xv(),amcisop->Yv(), amcisop->Zv());
1073 }
1074 else
1075 continue;
1076 }
1077 if(status!=1)
80c98845 1078 continue;
7a7744db 1079 if(mother == itrack)
c5f3236f 1080 continue;
7a7744db 1081 if(pt<ptcut)
80c98845 1082 continue;
80c98845 1083 if(vmcv.Perp()>1)
1084 continue;
7a7744db 1085 dR = TMath::Sqrt((phi-isophi)*(phi-isophi)+(eta-isoeta)*(eta-isoeta));
80c98845 1086 if(dR>radius)
1087 continue;
7a7744db 1088 sumpt += pt;
80c98845 1089 }
4b6ab175 1090 if(this->fDebug)
c5f3236f 1091 printf("\t\t::GetMcIsolation() returning value %f ...\n\n",sumpt);
80c98845 1092 return sumpt;
c5f3236f 1093 }
f1c66719 1094
33d8da67 1095//________________________________________________________________________
1096Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
1097{
1098 // Compute isolation based on tracks.
4b6ab175 1099 if(this->fDebug)
c5f3236f 1100 printf("\t::GetTrackIsolation() starting\n");
1101
33d8da67 1102 Double_t trkIsolation = 0;
1103 Double_t rad2 = radius*radius;
1104 Int_t ntrks = fSelPrimTracks->GetEntries();
1105 for(Int_t j = 0; j<ntrks; ++j) {
1106 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
1107 if (!track)
1108 continue;
1109 if (track->Pt()<pt)
1110 continue;
1111
1112 Float_t eta = track->Eta();
1113 Float_t phi = track->Phi();
1114 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1115 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1116 if(dist>rad2)
1117 continue;
1118 trkIsolation += track->Pt();
1119 }
4b6ab175 1120 if(this->fDebug)
c5f3236f 1121 printf("\t::GetTrackIsolation() returning\n\n");
33d8da67 1122 return trkIsolation;
1123}
f1c66719 1124
33d8da67 1125//________________________________________________________________________
1126Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const
1127{
f1c66719 1128 // Get phi band.
1129
33d8da67 1130 if(!fSelPrimTracks)
1131 return 0;
1132 Int_t ntracks = fSelPrimTracks->GetEntries();
1133 Double_t loweta = eta - R;
1134 Double_t upeta = eta + R;
1135 Double_t upphi = phi + R;
1136 Double_t lowphi = phi - R;
1137 Double_t et = 0;
1138 for(Int_t itr=0; itr<ntracks; itr++){
7a7744db 1139 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(itr));
33d8da67 1140 if(!track)
1141 continue;
1142 if(track->Pt()<minpt)
1143 continue;
1144 if((track->Eta() < upeta) && (track->Eta() > loweta))
1145 continue;
1146 if((track->Phi() > upphi) || (track->Phi() < lowphi))
1147 continue;
1148 et+=track->Pt();
1149 }
1150 return et;
1151}
f1c66719 1152
33d8da67 1153//________________________________________________________________________
1154Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
1155{
1156 // Calculate the energy of cross cells around the leading cell.
7a7744db 1157 if(!fVCells)
33d8da67 1158 return 0;
1159
33d8da67 1160 if (!fGeom)
1161 return 0;
1162
1163 Int_t iSupMod = -1;
1164 Int_t iTower = -1;
1165 Int_t iIphi = -1;
1166 Int_t iIeta = -1;
1167 Int_t iphi = -1;
1168 Int_t ieta = -1;
1169 Int_t iphis = -1;
1170 Int_t ietas = -1;
1171
1172 Double_t crossEnergy = 0;
1173
1174 fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
1175 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
1176
1177 Int_t ncells = cluster->GetNCells();
1178 for (Int_t i=0; i<ncells; i++) {
1179 Int_t cellAbsId = cluster->GetCellAbsId(i);
1180 fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1181 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1182 Int_t aphidiff = TMath::Abs(iphi-iphis);
1183 if (aphidiff>1)
1184 continue;
1185 Int_t aetadiff = TMath::Abs(ieta-ietas);
1186 if (aetadiff>1)
1187 continue;
1188 if ( (aphidiff==1 && aetadiff==0) ||
1189 (aphidiff==0 && aetadiff==1) ) {
7a7744db 1190 crossEnergy += fVCells->GetCellAmplitude(cellAbsId);
33d8da67 1191 }
1192 }
1193
1194 return crossEnergy;
1195}
1196
33d8da67 1197//________________________________________________________________________
1198Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1199{
1200 // Get maximum energy of attached cell.
1201
1202 id = -1;
7a7744db 1203 if(!fVCells)
33d8da67 1204 return 0;
1205
1206 Double_t maxe = 0;
1207 Int_t ncells = cluster->GetNCells();
1208 for (Int_t i=0; i<ncells; i++) {
7a7744db 1209 Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
33d8da67 1210 if (e>maxe) {
1211 maxe = e;
1212 id = cluster->GetCellAbsId(i);
1213 }
1214 }
1215 return maxe;
1216}
f1c66719 1217
33d8da67 1218//________________________________________________________________________
1219void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *)
1220{
33d8da67 1221 // Called once at the end of the query
ec681cf3 1222/* if(fIsGrid)
1223 return;*/
f1c66719 1224 if (fTree) {
c305511b 1225 printf("***tree %s being saved***\n",fTree->GetName());
55dcd0ae 1226 TFile *f = OpenFile(2);
33d8da67 1227 TDirectory::TContext context(f);
1228 if (f)
1229 fTree->Write();
1230 }
33d8da67 1231}