]> 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
changes from Prabhat Ranjan Pujahari <p.pujahari@cern.ch>
[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;
58acb8ab 670 if(fDebug)
671 printf("cluster %d survived\n", ic);
33d8da67 672 Float_t pos[3];
673 clus->GetPosition(pos);
674 TVector3 cpos(pos);
675 TString cellsAbsID;
b2d49404 676 AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++));
33d8da67 677 myclus->fE = clus->E();
678 myclus->fEt = clus->E()*TMath::Sin(cpos.Theta());
679 myclus->fR = cpos.Perp();
680 myclus->fEta = cpos.Eta();
681 myclus->fPhi = cpos.Phi();
682 if(cpos.Phi()<0){
683 myclus->fPhi+=TMath::TwoPi();
684 }
685 myclus->fN = clus->GetNCells();
686 Short_t id = -1;
687 myclus->fEmax = GetMaxCellEnergy( clus, id);
688 myclus->fIdmax = id;
7a7744db 689 myclus->fTmax = fVCells->GetCellTime(id);
33d8da67 690 myclus->fEcross = GetCrossEnergy( clus, id);
691 myclus->fDisp = clus->GetDispersion();
692 myclus->fM20 = clus->GetM20();
693 myclus->fM02 = clus->GetM02();
694 myclus->fTrDEta = clus->GetTrackDz();
695 myclus->fTrDPhi = clus->GetTrackDx();
696 myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
697 myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
698 myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
699 myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
700 myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.);
701 myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.);
702 myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.);
703 myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.);
704 for(Int_t icell=0;icell<myclus->fN;icell++){
705 Int_t absID = clus->GetCellAbsId(icell);
706 cellsAbsID.Append(Form("%d",absID));
707 cellsAbsID.Append(";");
708 }
709 myclus->fCellsAbsId = cellsAbsID;
710 myclus->fMcLabel = clus->GetLabel();
711 Int_t matchIndex = clus->GetTrackMatchedIndex();
7a7744db 712 if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){
33d8da67 713 myclus->fTrEp = -1;
714 continue;
715 }
7a7744db 716 AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex));
38faf770 717 if(!track)
33d8da67 718 continue;
7a7744db 719 if(fESD){
720 if(!fPrTrCuts)
721 continue;
722 if(!fPrTrCuts->IsSelected(track))
723 continue;
724 }
38faf770 725
33d8da67 726 myclus->fTrEp = clus->E()/track->P();
38faf770 727 myclus->fTrDedx = track->GetTPCsignal();
33d8da67 728 }
4b6ab175 729 if(this->fDebug)
c5f3236f 730 printf("::FillMyClusters() returning...\n\n");
33d8da67 731
732}
733//________________________________________________________________________
734void AliAnalysisTaskEMCALPhoton::FillMyAltClusters()
735{
f1c66719 736 // Fill clusters.
4b6ab175 737 if(this->fDebug)
c5f3236f 738 printf("::FillMyAltClusters() starting\n");
f1c66719 739
33d8da67 740 if(!fCaloClustersNew)
741 return;
742 Int_t nclus = fCaloClustersNew->GetEntries();
743 Int_t mcl = 0;
744 for(Int_t ic=0; ic < nclus; ic++){
7a7744db 745 AliVCluster *clus = static_cast<AliVCluster*>(fCaloClustersNew->At(ic));
33d8da67 746 if(!clus)
747 continue;
748 if(!clus->IsEMCAL())
749 continue;
750 if(clus->E() < fClusThresh)
751 continue;
752 Float_t pos[3];
753 clus->GetPosition(pos);
754 TVector3 cpos(pos);
755 TString cellsAbsID;
b2d49404 756 AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++));
33d8da67 757 myclus->fE = clus->E();
758 myclus->fEt = clus->E()*TMath::Sin(cpos.Theta());
759 myclus->fR = cpos.Perp();
760 myclus->fEta = cpos.Eta();
761 myclus->fPhi = cpos.Phi();
762 if(cpos.Phi()<0){
763 myclus->fPhi+=TMath::TwoPi();
764 }
765 myclus->fN = clus->GetNCells();
766 myclus->fDisp = clus->GetDispersion();
767 myclus->fM20 = clus->GetM20();
768 myclus->fM02 = clus->GetM02();
769 myclus->fTrDEta = clus->GetTrackDz();
770 myclus->fTrDPhi = clus->GetTrackDx();
771 myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
772 myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
773 myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
774 myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
775 for(Int_t icell=0;icell<myclus->fN;icell++){
776 Int_t absID = clus->GetCellAbsId(icell);
777 cellsAbsID.Append(Form("%d",absID));
778 cellsAbsID.Append(";");
779 }
780 myclus->fCellsAbsId = cellsAbsID;
781 myclus->fMcLabel = clus->GetLabel();
782 Int_t matchIndex = clus->GetTrackMatchedIndex();
7a7744db 783 if(matchIndex<0 || matchIndex>fVev->GetNumberOfTracks()){
33d8da67 784 myclus->fTrEp = -1;
785 continue;
786 }
7a7744db 787 AliVTrack* track = static_cast<AliVTrack*>(fVev->GetTrack(matchIndex));
33d8da67 788 if(!track){
789 myclus->fTrEp = -1;
790 continue;
791 }
792 if(!fPrTrCuts){
793 myclus->fTrEp = -1;
794 continue;
795 }
796 if(!fPrTrCuts->IsSelected(track)){
797 myclus->fTrEp = -1;
798 continue;
799 }
800 myclus->fTrEp = clus->E()/track->P();
801 }
4b6ab175 802 if(this->fDebug)
c5f3236f 803 printf("::FillMyAltClusters() returning...\n\n");
33d8da67 804
805}
806//________________________________________________________________________
3b37c011 807void AliAnalysisTaskEMCALPhoton::FillIsoTracks()
33d8da67 808{
f1c66719 809 // Fill high pt tracks.
4b6ab175 810 if(this->fDebug)
c5f3236f 811 printf("::FillIsoTracks() starting\n");
f1c66719 812
33d8da67 813 if(!fSelPrimTracks)
814 return;
815 Int_t ntracks = fSelPrimTracks->GetEntries();
816 Int_t imtr = 0;
817 for(Int_t it=0;it<ntracks; it++){
7a7744db 818 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(it));
33d8da67 819 if(!track)
820 continue;
01166399 821 /*if(track->Phi()<1.0 || track->Phi()>3.55)
822 continue;*/
3b37c011 823 if(TMath::Abs(track->Eta())>1.1)
824 continue;
a72860a2 825 /*if(track->Pt()<3)
826 continue;*/
827 AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
33d8da67 828 mtr->fPt = track->Pt();
829 mtr->fEta = track->Eta();
830 mtr->fPhi = track->Phi();
3b37c011 831 mtr->fCharge = track->Charge();
33d8da67 832 mtr->fDedx = track->GetTPCsignal();
833 mtr->fMcLabel = track->GetLabel();
834 }
4b6ab175 835 if(this->fDebug)
c5f3236f 836 printf("::FillIsoTracks() returning...\n\n");
33d8da67 837}
f1c66719 838
33d8da67 839//________________________________________________________________________
840void AliAnalysisTaskEMCALPhoton::GetMcParts()
841{
c5f3236f 842 // Get MC particles.
4b6ab175 843 if(this->fDebug)
c5f3236f 844 printf("::GetMcParts() starting\n");
f1c66719 845
4b6ab175 846 if (!this->fStack && !fAODMCParticles)
33d8da67 847 return;
4b6ab175 848 if(this->fDebug)
7a7744db 849 printf("either stack or aodmcpaticles exists\n");
850 const AliVVertex *evtVtx = 0;
4b6ab175 851 if(this->fStack)
7a7744db 852 evtVtx = fMCEvent->GetPrimaryVertex();
853 else
854 printf("no such thing as mc vvtx\n");
855 /*if (!evtVtx)
856 return;*/
4b6ab175 857 if(this->fDebug)
7a7744db 858 printf("mc vvertex ok\n");
859 Int_t nTracks = 0;
4b6ab175 860 if(this->fStack)
861 nTracks = this->fStack->GetNtrack();
7a7744db 862 else if(fAODMCParticles)
863 nTracks = fAODMCParticles->GetEntriesFast();
864 TParticle *mcP = 0;
865 AliAODMCParticle *amcP = 0;
4b6ab175 866 if(this->fDebug)
7a7744db 867 printf("loop in the %d mc particles starting\n",nTracks);
33d8da67 868 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
4b6ab175 869 if(this->fStack)
870 mcP = dynamic_cast<TParticle*>(this->fStack->Particle(iTrack));
7a7744db 871 if(fAODMCParticles)
872 amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack));
c5f3236f 873
33d8da67 874 // primary particle
7a7744db 875 Double_t dR = 0;
876 if(mcP){
877 dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
878 (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
879 (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
880 }
881 if((dR > 0.5))
33d8da67 882 continue;
c5f3236f 883
33d8da67 884
885 // kinematic cuts
7a7744db 886 Double_t pt = 0;
887 Double_t eta = 0;
888 Double_t phi = 0;
889 Int_t mother = -1;
890 Int_t pdgcode = 0;
891 if(mcP){
892 pt = mcP->Pt() ;
893 eta = mcP->Eta();
894 phi = mcP->Phi();
895 mother = mcP->GetMother(0);
896 pdgcode = mcP->GetPdgCode();
897 } else {
898 if(amcP){
899 pt = amcP->Pt();
900 eta = amcP->Eta();
901 phi = amcP->Phi();
902 mother = amcP->GetMother();
903 pdgcode = amcP->GetPdgCode();
904 } else
905 continue;
906 }
c5f3236f 907 if (pt<0.5)
33d8da67 908 continue;
c5f3236f 909
c5f3236f 910 if (TMath::Abs(eta)>0.7)
33d8da67 911 continue;
c5f3236f 912
c5f3236f 913 if (phi<1.0||phi>3.3)
33d8da67 914 continue;
c5f3236f 915
7a7744db 916 if (mother!=6 && mother!=7 )
c5f3236f 917 continue;
918
919
33d8da67 920 // pion or eta meson or direct photon
7a7744db 921 if(pdgcode == 111) {
922 } else if(pdgcode == 221) {
923 } else if(pdgcode == 22 ) {
80c98845 924 } else {
33d8da67 925 continue;
80c98845 926 }
c5f3236f 927
7a7744db 928 FillMcPart( fMyMcIndex++, iTrack);
33d8da67 929 }
4b6ab175 930 if(this->fDebug)
c5f3236f 931 printf("::GetMcParts() returning...\n\n");
33d8da67 932}
f1c66719 933
33d8da67 934//________________________________________________________________________
7a7744db 935void AliAnalysisTaskEMCALPhoton::FillMcPart( Int_t itrack, Int_t label)
33d8da67 936{
f1c66719 937 // Fill MC particles.
7a7744db 938 Int_t nTracks = 0;
939 TParticle *mcP = 0;
940 AliAODMCParticle *amcP= 0;
4b6ab175 941 if(this->fStack){
942 nTracks = this->fStack->GetNtrack();
943 mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack));
7a7744db 944 }
945 else if(fAODMCParticles){
946 nTracks = fAODMCParticles->GetEntriesFast();
947 amcP = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(itrack));
948 }
4b6ab175 949 if(this->fDebug)
c5f3236f 950 printf("\t::FillMcParts() starting with label %d\n", itrack);
7a7744db 951 TVector3 vmcv;
952 if(mcP)
953 vmcv.SetXYZ(mcP->Vx(),mcP->Vy(), mcP->Vz());
954 else {
955 if(amcP)
956 vmcv.SetXYZ(amcP->Xv(),amcP->Yv(), amcP->Zv());
957 else
958 return;
959 }
960
80c98845 961 AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(itrack));
c5f3236f 962 mcp->fLabel = label ;
7a7744db 963 if(mcP){
964 mcp->fPdg = mcP->GetPdgCode() ;
965 mcp->fPt = mcP->Pt() ;
966 mcp->fEta = mcP->Eta() ;
967 mcp->fPhi = mcP->Phi() ;
968 mcp->fMother = mcP->GetMother(0) ;
969 mcp->fFirstD = mcP->GetFirstDaughter() ;
970 mcp->fLastD = mcP->GetLastDaughter() ;
971 mcp->fStatus = mcP->GetStatusCode();
972 }
973 if(amcP){
974 mcp->fPdg = amcP->GetPdgCode() ;
975 mcp->fPt = amcP->Pt() ;
976 mcp->fEta = amcP->Eta() ;
977 mcp->fPhi = amcP->Phi() ;
978 mcp->fMother = amcP->GetMother() ;
979 mcp->fFirstD = amcP->GetDaughter(0) ;
980 mcp->fLastD = amcP->GetDaughter(amcP->GetNDaughters()-1) ;
981 mcp->fStatus = amcP->GetStatus();
982 }
33d8da67 983 mcp->fVR = vmcv.Perp();
984 if(vmcv.Perp()>0){
985 mcp->fVEta = vmcv.Eta() ;
986 mcp->fVPhi = vmcv.Phi() ;
987 }
7a7744db 988 if(itrack == 8){
989 mcp->fIso = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.4 , 0.2);
990 mcp->fIso3 = AliAnalysisTaskEMCALPhoton::GetMcIsolation( itrack, 0.3 , 0.2);
991 }
4b6ab175 992 if(this->fDebug)
7a7744db 993 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);
994 for(Int_t id=mcp->fFirstD; id < mcp->fLastD; id++){
995 if(id<=mcp->fMother)
c5f3236f 996 continue;
997 if(id<0 || id>nTracks)
998 continue;
7a7744db 999 FillMcPart( fMyMcIndex++, id);
c5f3236f 1000 }
1001
80c98845 1002}
1003//________________________________________________________________________
7a7744db 1004Double_t AliAnalysisTaskEMCALPhoton::GetMcIsolation( Int_t itrack, Double_t radius, Double_t ptcut) const
80c98845 1005{
4b6ab175 1006 if(this->fDebug){
c5f3236f 1007 printf("\t\t::GetMcIsolation() starting\n");
7a7744db 1008 //printf("\t\t incoming particle: PDG = %d, itrack=%d;\n",mcP->GetPdgCode(),itrack);
c5f3236f 1009 }
4b6ab175 1010 if (!this->fStack && !this->fAODMCParticles){
7a7744db 1011 printf("\t\t\tNo MC stack/array!\n");
80c98845 1012 return -1;
7a7744db 1013 }
80c98845 1014 if(itrack<6 || itrack>8)
1015 return -1;
4b6ab175 1016 if(this->fDebug)
7a7744db 1017 printf("\t\t\tparticle of interest selected\n");
1018 TParticle *mcP = 0;
1019 AliAODMCParticle *amcP = 0;
1020 Int_t pdgcode = 0;
1021 Float_t eta = 0;
1022 Float_t phi = 0;
80c98845 1023 Double_t sumpt=0;
80c98845 1024 Float_t dR;
7a7744db 1025 Int_t nparts = 0;
4b6ab175 1026 if(this->fStack){
7a7744db 1027 nparts = fStack->GetNtrack();
4b6ab175 1028 mcP = dynamic_cast<TParticle*>(this->fStack->Particle(itrack));
7a7744db 1029 eta = mcP->Eta();
1030 phi = mcP->Phi();
1031 pdgcode = mcP->GetPdgCode();
1032 }
4b6ab175 1033 if(this->fAODMCParticles){
7a7744db 1034 nparts = fAODMCParticles->GetEntriesFast();
4b6ab175 1035 amcP = dynamic_cast<AliAODMCParticle*>(this->fAODMCParticles->At(itrack));
1036 if(amcP){
1037 eta = amcP->Eta();
1038 phi = amcP->Phi();
1039 pdgcode = amcP->GetPdgCode();
1040 }
7a7744db 1041 }
1042 if(pdgcode!=22)
1043 return -1;
1044 TParticle *mcisop = 0;
1045 AliAODMCParticle *amcisop = 0;
80c98845 1046 for(Int_t ip = 0; ip<nparts; ip++){
80c98845 1047 if(ip==itrack)
1048 continue;
4b6ab175 1049 if(this->fStack)
1050 mcisop = dynamic_cast<TParticle*>(this->fStack->Particle(ip));
7a7744db 1051 if(fAODMCParticles)
1052 amcisop = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(ip));
1053 Int_t status = 0;
1054 Int_t mother = 0;
1055 Float_t pt = 0;
1056 Float_t isophi = -99;
1057 Float_t isoeta = -99;
1058 TVector3 vmcv;
1059 if(mcisop){
1060 status = mcisop->GetStatusCode();
1061 mother = mcisop->GetMother(0);
1062 pt = mcisop->Pt();
1063 isophi = mcisop->Phi();
1064 isoeta = mcisop->Eta();
1065 vmcv.SetXYZ(mcisop->Vx(),mcisop->Vy(), mcisop->Vz());
1066 }
1067 else {
1068 if(amcisop){
1069 status = amcisop->GetStatus();
1070 mother = amcisop->GetMother();
1071 pt = amcisop->Pt();
1072 isophi = amcisop->Phi();
1073 isoeta = amcisop->Eta();
1074 vmcv.SetXYZ(amcisop->Xv(),amcisop->Yv(), amcisop->Zv());
1075 }
1076 else
1077 continue;
1078 }
1079 if(status!=1)
80c98845 1080 continue;
7a7744db 1081 if(mother == itrack)
c5f3236f 1082 continue;
7a7744db 1083 if(pt<ptcut)
80c98845 1084 continue;
80c98845 1085 if(vmcv.Perp()>1)
1086 continue;
7a7744db 1087 dR = TMath::Sqrt((phi-isophi)*(phi-isophi)+(eta-isoeta)*(eta-isoeta));
80c98845 1088 if(dR>radius)
1089 continue;
7a7744db 1090 sumpt += pt;
80c98845 1091 }
4b6ab175 1092 if(this->fDebug)
c5f3236f 1093 printf("\t\t::GetMcIsolation() returning value %f ...\n\n",sumpt);
80c98845 1094 return sumpt;
c5f3236f 1095 }
f1c66719 1096
33d8da67 1097//________________________________________________________________________
1098Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
1099{
1100 // Compute isolation based on tracks.
4b6ab175 1101 if(this->fDebug)
c5f3236f 1102 printf("\t::GetTrackIsolation() starting\n");
1103
33d8da67 1104 Double_t trkIsolation = 0;
1105 Double_t rad2 = radius*radius;
1106 Int_t ntrks = fSelPrimTracks->GetEntries();
1107 for(Int_t j = 0; j<ntrks; ++j) {
1108 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
1109 if (!track)
1110 continue;
1111 if (track->Pt()<pt)
1112 continue;
1113
1114 Float_t eta = track->Eta();
1115 Float_t phi = track->Phi();
1116 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1117 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1118 if(dist>rad2)
1119 continue;
1120 trkIsolation += track->Pt();
1121 }
4b6ab175 1122 if(this->fDebug)
c5f3236f 1123 printf("\t::GetTrackIsolation() returning\n\n");
33d8da67 1124 return trkIsolation;
1125}
f1c66719 1126
33d8da67 1127//________________________________________________________________________
1128Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const
1129{
f1c66719 1130 // Get phi band.
1131
33d8da67 1132 if(!fSelPrimTracks)
1133 return 0;
1134 Int_t ntracks = fSelPrimTracks->GetEntries();
1135 Double_t loweta = eta - R;
1136 Double_t upeta = eta + R;
1137 Double_t upphi = phi + R;
1138 Double_t lowphi = phi - R;
1139 Double_t et = 0;
1140 for(Int_t itr=0; itr<ntracks; itr++){
7a7744db 1141 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(itr));
33d8da67 1142 if(!track)
1143 continue;
1144 if(track->Pt()<minpt)
1145 continue;
1146 if((track->Eta() < upeta) && (track->Eta() > loweta))
1147 continue;
1148 if((track->Phi() > upphi) || (track->Phi() < lowphi))
1149 continue;
1150 et+=track->Pt();
1151 }
1152 return et;
1153}
f1c66719 1154
33d8da67 1155//________________________________________________________________________
1156Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
1157{
1158 // Calculate the energy of cross cells around the leading cell.
7a7744db 1159 if(!fVCells)
33d8da67 1160 return 0;
1161
33d8da67 1162 if (!fGeom)
1163 return 0;
1164
1165 Int_t iSupMod = -1;
1166 Int_t iTower = -1;
1167 Int_t iIphi = -1;
1168 Int_t iIeta = -1;
1169 Int_t iphi = -1;
1170 Int_t ieta = -1;
1171 Int_t iphis = -1;
1172 Int_t ietas = -1;
1173
1174 Double_t crossEnergy = 0;
1175
1176 fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
1177 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
1178
1179 Int_t ncells = cluster->GetNCells();
1180 for (Int_t i=0; i<ncells; i++) {
1181 Int_t cellAbsId = cluster->GetCellAbsId(i);
1182 fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
1183 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
1184 Int_t aphidiff = TMath::Abs(iphi-iphis);
1185 if (aphidiff>1)
1186 continue;
1187 Int_t aetadiff = TMath::Abs(ieta-ietas);
1188 if (aetadiff>1)
1189 continue;
1190 if ( (aphidiff==1 && aetadiff==0) ||
1191 (aphidiff==0 && aetadiff==1) ) {
7a7744db 1192 crossEnergy += fVCells->GetCellAmplitude(cellAbsId);
33d8da67 1193 }
1194 }
1195
1196 return crossEnergy;
1197}
1198
33d8da67 1199//________________________________________________________________________
1200Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1201{
1202 // Get maximum energy of attached cell.
1203
1204 id = -1;
7a7744db 1205 if(!fVCells)
33d8da67 1206 return 0;
1207
1208 Double_t maxe = 0;
1209 Int_t ncells = cluster->GetNCells();
1210 for (Int_t i=0; i<ncells; i++) {
7a7744db 1211 Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
33d8da67 1212 if (e>maxe) {
1213 maxe = e;
1214 id = cluster->GetCellAbsId(i);
1215 }
1216 }
1217 return maxe;
1218}
f1c66719 1219
33d8da67 1220//________________________________________________________________________
1221void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *)
1222{
33d8da67 1223 // Called once at the end of the query
ec681cf3 1224/* if(fIsGrid)
1225 return;*/
f1c66719 1226 if (fTree) {
c305511b 1227 printf("***tree %s being saved***\n",fTree->GetName());
55dcd0ae 1228 TFile *f = OpenFile(2);
33d8da67 1229 TDirectory::TContext context(f);
1230 if (f)
1231 fTree->Write();
1232 }
33d8da67 1233}