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