]> 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
make work with AOD
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALPhoton.cxx
CommitLineData
33d8da67 1#include "TChain.h"
2#include "TTree.h"
3#include "TH1F.h"
4#include "TH2F.h"
5#include "TCanvas.h"
6
7#include "AliAnalysisTask.h"
8#include "AliAnalysisManager.h"
9
10#include "AliESDEvent.h"
11#include "AliESDHeader.h"
12#include "AliESDUtils.h"
13#include "AliESDInputHandler.h"
14#include "AliCentrality.h"
15#include "AliESDpid.h"
16#include "AliKFParticle.h"
17
18#include "AliMCEventHandler.h"
19#include "AliMCEvent.h"
20#include "AliStack.h"
21#include "TParticle.h"
22
23
24#include "AliESDtrackCuts.h"
25#include "AliESDv0.h"
26#include "AliV0vertexer.h"
27#include "AliESDCaloCluster.h"
28#include "AliESDCaloCells.h"
29#include "AliEMCALGeometry.h"
30#include "AliEMCALRecoUtils.h"
31#include "AliVCluster.h"
32#include "AliAnalysisTaskEMCALClusterizeFast.h"
33#include "TLorentzVector.h"
34
35
36#include "AliAnalysisTaskEMCALPhoton.h"
37#include "TFile.h"
38
39
40ClassImp(AliAnalysisTaskEMCALPhoton)
41
42//________________________________________________________________________
43AliAnalysisTaskEMCALPhoton::AliAnalysisTaskEMCALPhoton(const char *name)
44 : AliAnalysisTaskSE(name),
45 fTrCuts(0),
46 fPrTrCuts(0),
47 fSelTracks(0),
48 fSelPrimTracks(0),
49 fPhotConvArray(0),
50 fMyClusts(0),
51 fMyAltClusts(0),
52 fMyCells(0),
53 fMyTracks(0),
54 fMyMcParts(0),
55 fHeader(0),
56 fCaloClusters(0),
57 fCaloClustersNew(0),
58 fEMCalCells(0),
59 fGeom(0x0),
60 fTimeResTOF(0),
61 fMipResponseTPC(0),
62 fGeoName("EMCAL_COMPLETEV1"),
63 fPeriod("LHC11c"),
64 fIsTrain(0),
65 fIsMC(0),
66 fIsGrid(0),
67 fClusThresh(2.),
68 fClusterizer(0),
69 fCaloClustersName("EmcalClusterv2"),
70
71
72 fESD(0),
73 fMCEvent(0),
74 fStack(0x0),
75
76 fOutputList(0),
77 fTree(0),
78
79
80 fNV0sBefAndAftRerun(0),
81 fConversionVtxXY(0),
82 fInvMassV0(0),
83 fInvMassV0KF(0),
84 fInvMassV0SS(0),
85 fDedxPAll(0)
86
87{
88 // Constructor
89
90 // Define input and output slots here
91 // Input slot #0 works with a TChain
92 DefineInput(0, TChain::Class());
93 // Output slot #0 id reserved by the base class for AOD
94 // Output slot #1 writes into a TH1 container
95 DefineOutput(1, TList::Class());
96}
97//________________________________________________________________________
98void AliAnalysisTaskEMCALPhoton::UserCreateOutputObjects()
99{
100 // Create histograms
101 // Called once
102
103 fSelTracks = new TObjArray();
104
105 fSelPrimTracks = new TObjArray();
106
b2d49404 107 if (TClass::GetClass("AliPhotonHeaderObj"))
108 TClass::GetClass("AliPhotonHeaderObj")->IgnoreTObjectStreamer();
109 fHeader = new AliPhotonHeaderObj;
33d8da67 110
b2d49404 111 if (TClass::GetClass("AliPhotonConvObj"))
112 TClass::GetClass("AliPhotonConvObj")->IgnoreTObjectStreamer();
113 fPhotConvArray = new TClonesArray("AliPhotonConvObj");
33d8da67 114
b2d49404 115 if (TClass::GetClass("AliPhotonClusterObj"))
116 TClass::GetClass("AliPhotonClusterObj")->IgnoreTObjectStreamer();
117 fMyClusts = new TClonesArray("AliPhotonClusterObj");
33d8da67 118
b2d49404 119 if (TClass::GetClass("AliPhotonCellObj"))
120 TClass::GetClass("AliPhotonCellObj")->IgnoreTObjectStreamer();
121 fMyCells = new TClonesArray("AliPhotonCellObj");
33d8da67 122
b2d49404 123 if (TClass::GetClass("AliPhotonTrackObj"))
124 TClass::GetClass("AliPhotonTrackObj")->IgnoreTObjectStreamer();
125 fMyTracks = new TClonesArray("AliPhotonTrackObj");
33d8da67 126
127 if (fClusterizer || fIsTrain){
128 if(fClusterizer)
129 fCaloClustersName = fClusterizer->GetNewClusterArrayName();
130 else {
131 if(fPeriod.Contains("10h") || fPeriod.Contains("11h"))
132 fCaloClustersName = "EmcalClustersv1";
133 else
134 fCaloClustersName = "EmcalClustersv2";
135 }
b2d49404 136 if (TClass::GetClass("AliPhotonClusterObj"))
137 TClass::GetClass("AliPhotonClusterObj")->IgnoreTObjectStreamer();
138 fMyAltClusts = new TClonesArray("AliPhotonClusterObj");
33d8da67 139 }
140 cout<<fCaloClustersName<<endl;
141 if(fIsMC){
b2d49404 142 if (TClass::GetClass("AliPhotonMcPartObj"))
143 TClass::GetClass("AliPhotonMcPartObj")->IgnoreTObjectStreamer();
144 fMyMcParts = new TClonesArray("AliPhotonMcPartObj");
33d8da67 145 }
146
147 fCaloClusters = new TRefArray();
148
149 fOutputList = new TList();
150 fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging)
151
152 if( !fTree){
153 TFile *f = OpenFile(1);
154 TDirectory::TContext context(f);
155 if (f) {
156 f->SetCompressionLevel(2);
157 fTree = new TTree("photon_ana_out", "StandaloneTree");
158 fTree->SetDirectory(f);
159 if (fIsTrain) {
160 fTree->SetAutoFlush(-2*1024*1024);
161 fTree->SetAutoSave(0);
162 } else {
163 fTree->SetAutoFlush(-32*1024*1024);
164 fTree->SetAutoSave(0);
165 }
166
167 fTree->Branch("header", &fHeader, 16*1024, 99);
168 fTree->Branch("conversions", &fPhotConvArray, 8*16*1024, 99);
169 fTree->Branch("clusters", &fMyClusts, 8*16*1024, 99);
170 if(fClusterizer || fIsTrain)
171 fTree->Branch(fCaloClustersName, &fMyAltClusts, 8*16*1024, 99);
172 fTree->Branch("cells", &fMyCells, 8*16*1024, 99);
173 fTree->Branch("HighPtTracks", &fMyTracks, 8*16*1024, 99);
174 if(fIsMC)
175 fTree->Branch("mcparts", &fMyMcParts, 8*16*1024, 99);
176 }
177 }
178 if(fIsGrid)fOutputList->Add(fTree);
179 fGeom = AliEMCALGeometry::GetInstance(fGeoName);
180
181
182 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);
183 fOutputList->Add(fNV0sBefAndAftRerun);
184
185 fConversionVtxXY = new TH2F("hConversionVtxXY","x and y of conversion vertex candidates;x;y",1000,-100,100,1000,-100,100);
186 fOutputList->Add(fConversionVtxXY);
187
188 fInvMassV0 = new TH1F("hInvMassV0","v0->GetEffMass();v0->GetEffMass();dN/dM",400,0,4);
189 fOutputList->Add(fInvMassV0);
190
191 fInvMassV0KF = new TH1F("hInvMassV0KF","Inv. mass calculated from AliKFParticle made of V0 tracks;mass_{TrTr};dN/dM",400,0,4);
192 fOutputList->Add(fInvMassV0KF);
193
194 fInvMassV0SS = new TH1F("hInvMassV0SS","Inv. mass (same sign) calculated from AliKFParticle made of V0 tracks;mass_{TrTr};dN/dM",400,0,4);
195 fOutputList->Add(fInvMassV0SS);
196
197 fDedxPAll = new TH2F("hDedxPAll","dE/dx vs p (all selected tracks);p (GeV/c);dE/dx (a.u.)",400,0,40, 150,0,150);
198 fOutputList->Add(fDedxPAll);
199
200
201 PostData(1, fOutputList);
202}
203
204//________________________________________________________________________
205void AliAnalysisTaskEMCALPhoton::UserExec(Option_t *)
206{
207 //event trigger selection
208 Bool_t isSelected = 0;
209 if(fPeriod.Contains("11")){
210 if(fPeriod.Contains("11a"))
211 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1);
212 if(fPeriod.Contains("11c") ||fPeriod.Contains("11d") )
213 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7);
214 if(fPeriod.Contains("11h") )
215 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA);
216
217 }
218 else
219 isSelected = 1;
220/* if(!isSelected )
221 return; */
222 // Main loop
223 // Called for each event
224
225 // Post output data.
226 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
227 if (!fESD) {
228 printf("ERROR: fESD not available, returning...\n");
229 return;
230 }
231
232
233
234 AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
235 if(!pv)
236 return;
237 if(TMath::Abs(pv->GetZ())>15)
238 return;
239 //fCaloClustersNew = dynamic_cast<TClonesArray*>(fESD->FindListObject(fCaloClustersName));
240
241
242 // Track loop to fill a pT spectrum
243 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
244 AliESDtrack* track = fESD->GetTrack(iTracks);
245 if (!track)
246 continue;
247
248
249 if (fTrCuts && fTrCuts->IsSelected(track)){
250 fSelTracks->Add(track);
251 fDedxPAll->Fill(track->P(), track->GetTPCsignal());
252 }
253 if (fPrTrCuts && fPrTrCuts->IsSelected(track))
254 fSelPrimTracks->Add(track);
255
256
257
258 } //track loop
259 fHeader->fTrClassMask = fESD->GetHeader()->GetTriggerMask();
260 fHeader->fTrCluster = fESD->GetHeader()->GetTriggerCluster();
261 AliCentrality *cent = InputEvent()->GetCentrality();
262 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
263 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
264 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
265 if(!fIsTrain){
266 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
267 if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
268 break;
269 fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
270 }
271 }
272 fESD->GetEMCALClusters(fCaloClusters);
273 fHeader->fNClus = fCaloClusters->GetEntries();
274 fEMCalCells = fESD->GetEMCALCells();
275 fHeader->fNCells = fEMCalCells->GetNumberOfCells();
276 fMCEvent = MCEvent();
277 if(fMCEvent)
278 fStack = (AliStack*)fMCEvent->Stack();
279
280
281 FindConversions();
282// cout<<"****calling FillMyCells()...\n";
283 FillMyCells();
284// cout<<"****calling FillMyClus...\n";
285 FillMyClusters();
286 if(fCaloClustersNew)
287 FillMyAltClusters();
288// cout<<"****calling FillHighPt...\n";
289 FillHighPtTracks();
290 if(fIsMC)
291 GetMcParts();
292
293 fTree->Fill();
294 fSelTracks->Clear();
295 fSelPrimTracks->Clear();
296 fPhotConvArray->Clear();
297 fMyClusts->Clear();
298 if(fMyAltClusts)
299 fMyAltClusts->Clear();
300 fMyCells->Clear();
301 fMyTracks->Clear();
302 if(fMyMcParts)
303 fMyMcParts->Clear();
304 fCaloClusters->Clear();
305 if(fCaloClustersNew)
306 fCaloClustersNew->Clear();
307 PostData(1, fOutputList);
308}
309
310//________________________________________________________________________
311void AliAnalysisTaskEMCALPhoton::FindConversions()
312{
313 if(!fTrCuts)
314 return;
315 Int_t iconv = 0;
316 Int_t nV0Orig = fESD->GetNumberOfV0s();
317 Int_t ntracks = fESD->GetNumberOfTracks();
318 fESD->ResetV0s();
319 AliV0vertexer lV0vtxer;
320 lV0vtxer.Tracks2V0vertices(fESD);
321 Int_t nV0s = fESD->GetNumberOfV0s();
322 fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s);
323 for(Int_t iv0=0; iv0<nV0s; iv0++){
324 AliESDv0 * v0 = fESD->GetV0(iv0);
325 if(!v0)
326 continue;
327 Double_t vpos[3];
328 fInvMassV0->Fill(v0->GetEffMass());
329 v0->GetXYZ(vpos[0], vpos[1], vpos[2]);
330 Int_t ipos = v0->GetPindex();
331 Int_t ineg = v0->GetNindex();
332 if(ipos<0 || ipos> ntracks)
333 continue;
334 if(ineg<0 || ineg> ntracks)
335 continue;
336 AliESDtrack *pos = static_cast<AliESDtrack*>(fESD->GetTrack(ipos));
337 AliESDtrack *neg = static_cast<AliESDtrack*>(fESD->GetTrack(ineg));
338 if(!pos || !neg)
339 continue;
340 if (!fTrCuts->IsSelected(pos) || !fTrCuts->IsSelected(neg))
341 continue;
342 /*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65)
343 continue;*/
344 const AliExternalTrackParam * paramPos = v0->GetParamP() ;
345 const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
346 if(!paramPos || !paramNeg)
347 continue;
348 if(pos->GetSign() <0){//change tracks
349 pos=neg ;
350 neg=fESD->GetTrack(v0->GetPindex()) ;
351 paramPos=paramNeg ;
352 paramNeg=v0->GetParamP() ;
353 }
354 AliKFParticle negKF(*paramNeg,11);
355 AliKFParticle posKF(*paramPos,-11);
356 AliKFParticle photKF(negKF,posKF) ;
357
358
359 if( neg->GetKinkIndex(0) > 0 || pos->GetKinkIndex(0) > 0)
360 continue ;
361
362
363 if(pos->GetSign() == neg->GetSign()){
364 fInvMassV0SS->Fill(photKF.GetMass());
365 continue;
366 }
367 fInvMassV0KF->Fill(photKF.GetMass());
368 TLorentzVector photLV(photKF.GetPx(), photKF.GetPy(),photKF.GetPz(), photKF.GetE());
369 if(photLV.M()>0.05 || photLV.M()<0)
370 continue;
371 fConversionVtxXY->Fill(vpos[0], vpos[1]);//negKF.GetX(), negKF.GetY());
372 Double_t convPhi = TMath::ATan2(photLV.Py(),photLV.Px());
373 if(convPhi<0)
374 convPhi+=TMath::TwoPi();
375 TVector3 vecpos(vpos);
376 Double_t v0Phi = 0;
377 if(vecpos.Perp()>0)
378 vecpos.Phi();
379 if(v0Phi<0)
380 v0Phi+=TMath::TwoPi();
b2d49404 381 AliPhotonConvObj *myconv = static_cast<AliPhotonConvObj*>(fPhotConvArray->New(iconv++));
33d8da67 382 myconv->fPt = photLV.Pt();
383 myconv->fEta = photLV.Eta();
384 myconv->fPhi = convPhi;
385 myconv->fVR = vecpos.Perp();
386 if(vecpos.Perp()>0)
387 myconv->fVEta = vecpos.Eta();
388 myconv->fVPhi = v0Phi;
389 myconv->fMass = photLV.M();
390 myconv->fMcLabel = -3; //WARNING: include the correct labeling
391 //negative daughter
392 myconv->fNegPt = negKF.GetPt();
393 myconv->fNegEta = negKF.GetEta();
394 Double_t trackPhi = negKF.GetPhi();
395 if(trackPhi<0)
396 trackPhi+=TMath::TwoPi();
397 myconv->fNegPhi = trackPhi;
398 myconv->fNegDedx = neg->GetTPCsignal();
399 myconv->fNegMcLabel = neg->GetLabel();
400 //negative daughter
401 myconv->fPosPt = posKF.GetPt();
402 myconv->fPosEta = posKF.GetEta();
403 trackPhi = posKF.GetPhi();
404 if(trackPhi<0)
405 trackPhi+=TMath::TwoPi();
406 myconv->fPosPhi = trackPhi;
407 myconv->fPosDedx = pos->GetTPCsignal();
408 myconv->fPosMcLabel = pos->GetLabel();
409 }
410 return;
411}
412//________________________________________________________________________
413void AliAnalysisTaskEMCALPhoton::FillMyCells()
414{
415 if(!fEMCalCells)
416 return;
417// cout<<"++++FillMyCells starts!\n";
418 Int_t ncells = fEMCalCells->GetNumberOfCells();
419 Int_t mcel = 0;
420 for(Int_t icell = 0; icell<ncells; icell++){
421 Int_t absID = TMath::Abs(fEMCalCells->GetCellNumber(icell));
b2d49404 422 AliPhotonCellObj *mycell = static_cast<AliPhotonCellObj*>(fMyCells->New(mcel++));
33d8da67 423 Float_t eta=-1, phi=-1;
424 if(!fGeom){
425 cout<<"+++fGeom not found! MyCells branch will not be filled for this event!+++\n";
426 return;
427 }
428 if(!fGeom)
429 return;
430 if(!fIsMC)fGeom->EtaPhiFromIndex(absID,eta,phi);
431 Float_t theta = 2*TMath::ATan(TMath::Exp(-eta));
432 mycell->fAbsID = absID;
433 mycell->fE = fEMCalCells->GetCellAmplitude(absID);
434 mycell->fEt = fEMCalCells->GetCellAmplitude(absID)*TMath::Sin(theta);
435 mycell->fEta = eta;
436 mycell->fPhi = phi;
437 mycell->fTime = fEMCalCells->GetCellTime(absID);
438 }
439// cout<<"++++FillMyCells done!\n";
440}
441//________________________________________________________________________
442void AliAnalysisTaskEMCALPhoton::FillMyClusters()
443{
444 if(!fCaloClusters)
445 return;
446 Int_t nclus = fCaloClusters->GetEntries();
447 Int_t mcl = 0;
448 for(Int_t ic=0; ic < nclus; ic++){
449 AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
450 if(!clus)
451 continue;
452 if(!clus->IsEMCAL())
453 continue;
454 if(clus->E() < fClusThresh)
455 continue;
456 Float_t pos[3];
457 clus->GetPosition(pos);
458 TVector3 cpos(pos);
459 TString cellsAbsID;
b2d49404 460 AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++));
33d8da67 461 myclus->fE = clus->E();
462 myclus->fEt = clus->E()*TMath::Sin(cpos.Theta());
463 myclus->fR = cpos.Perp();
464 myclus->fEta = cpos.Eta();
465 myclus->fPhi = cpos.Phi();
466 if(cpos.Phi()<0){
467 myclus->fPhi+=TMath::TwoPi();
468 }
469 myclus->fN = clus->GetNCells();
470 Short_t id = -1;
471 myclus->fEmax = GetMaxCellEnergy( clus, id);
472 myclus->fIdmax = id;
473 myclus->fTmax = fEMCalCells->GetCellTime(id);
474 myclus->fEcross = GetCrossEnergy( clus, id);
475 myclus->fDisp = clus->GetDispersion();
476 myclus->fM20 = clus->GetM20();
477 myclus->fM02 = clus->GetM02();
478 myclus->fTrDEta = clus->GetTrackDz();
479 myclus->fTrDPhi = clus->GetTrackDx();
480 myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
481 myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
482 myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
483 myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
484 myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.);
485 myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.);
486 myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.);
487 myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.);
488 for(Int_t icell=0;icell<myclus->fN;icell++){
489 Int_t absID = clus->GetCellAbsId(icell);
490 cellsAbsID.Append(Form("%d",absID));
491 cellsAbsID.Append(";");
492 }
493 myclus->fCellsAbsId = cellsAbsID;
494 myclus->fMcLabel = clus->GetLabel();
495 Int_t matchIndex = clus->GetTrackMatchedIndex();
496 if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
497 myclus->fTrEp = -1;
498 continue;
499 }
500 AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
501 if(!track){
502 myclus->fTrEp = -1;
503 continue;
504 }
505 if(!fPrTrCuts){
506 myclus->fTrEp = -1;
507 continue;
508 }
509 if(!fPrTrCuts->IsSelected(track)){
510 myclus->fTrEp = -1;
511 continue;
512 }
513 myclus->fTrEp = clus->E()/track->P();
514 }
515
516}
517//________________________________________________________________________
518void AliAnalysisTaskEMCALPhoton::FillMyAltClusters()
519{
520 if(!fCaloClustersNew)
521 return;
522 Int_t nclus = fCaloClustersNew->GetEntries();
523 Int_t mcl = 0;
524 for(Int_t ic=0; ic < nclus; ic++){
525 AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClustersNew->At(ic));
526 if(!clus)
527 continue;
528 if(!clus->IsEMCAL())
529 continue;
530 if(clus->E() < fClusThresh)
531 continue;
532 Float_t pos[3];
533 clus->GetPosition(pos);
534 TVector3 cpos(pos);
535 TString cellsAbsID;
b2d49404 536 AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++));
33d8da67 537 myclus->fE = clus->E();
538 myclus->fEt = clus->E()*TMath::Sin(cpos.Theta());
539 myclus->fR = cpos.Perp();
540 myclus->fEta = cpos.Eta();
541 myclus->fPhi = cpos.Phi();
542 if(cpos.Phi()<0){
543 myclus->fPhi+=TMath::TwoPi();
544 }
545 myclus->fN = clus->GetNCells();
546 myclus->fDisp = clus->GetDispersion();
547 myclus->fM20 = clus->GetM20();
548 myclus->fM02 = clus->GetM02();
549 myclus->fTrDEta = clus->GetTrackDz();
550 myclus->fTrDPhi = clus->GetTrackDx();
551 myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
552 myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
553 myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
554 myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
555 for(Int_t icell=0;icell<myclus->fN;icell++){
556 Int_t absID = clus->GetCellAbsId(icell);
557 cellsAbsID.Append(Form("%d",absID));
558 cellsAbsID.Append(";");
559 }
560 myclus->fCellsAbsId = cellsAbsID;
561 myclus->fMcLabel = clus->GetLabel();
562 Int_t matchIndex = clus->GetTrackMatchedIndex();
563 if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
564 myclus->fTrEp = -1;
565 continue;
566 }
567 AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
568 if(!track){
569 myclus->fTrEp = -1;
570 continue;
571 }
572 if(!fPrTrCuts){
573 myclus->fTrEp = -1;
574 continue;
575 }
576 if(!fPrTrCuts->IsSelected(track)){
577 myclus->fTrEp = -1;
578 continue;
579 }
580 myclus->fTrEp = clus->E()/track->P();
581 }
582
583}
584//________________________________________________________________________
585void AliAnalysisTaskEMCALPhoton::FillHighPtTracks()
586{
587 if(!fSelPrimTracks)
588 return;
589 Int_t ntracks = fSelPrimTracks->GetEntries();
590 Int_t imtr = 0;
591 for(Int_t it=0;it<ntracks; it++){
592 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(it));
593 if(!track)
594 continue;
595 if(track->Pt()<3)
596 continue;
b2d49404 597 AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
33d8da67 598 mtr->fPt = track->Pt();
599 mtr->fEta = track->Eta();
600 mtr->fPhi = track->Phi();
601 mtr->fDedx = track->GetTPCsignal();
602 mtr->fMcLabel = track->GetLabel();
603 }
604}
605//________________________________________________________________________
606void AliAnalysisTaskEMCALPhoton::GetMcParts()
607{
608 if (!fStack)
609 return;
610// cout<<"++++Get - MC starts!\n";
611
612 const AliVVertex *evtVtx = fMCEvent->GetPrimaryVertex();
613 if (!evtVtx)
614 return;
615 Int_t ipart = 0;
616 Int_t nTracks = fStack->GetNtrack();
617 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
618 TParticle *mcP = static_cast<TParticle*>(fStack->Particle(iTrack));
619 if (!mcP)
620 continue;
621 // primary particle
622 Double_t dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
623 (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
624 (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
625 if(dR > 0.5)
626 continue;
627
628 // kinematic cuts
629 Double_t pt = mcP->Pt() ;
630 if (pt<0.5)
631 continue;
632 Double_t eta = mcP->Eta();
633 if (eta<-0.7||eta>0.7)
634 continue;
635 Double_t phi = mcP->Phi();
636 if (phi<1.0||phi>3.3)
637 continue;
638 // pion or eta meson or direct photon
639 if(mcP->GetPdgCode() == 111) {
640 } else if(mcP->GetPdgCode() == 221) {
641 } else if(mcP->GetPdgCode() == 22 ) {
642 } else
643 continue;
644 bool checkIfAlreadySaved = false;
645 for(Int_t imy=0;imy<fMyMcParts->GetEntries();imy++){
b2d49404 646 AliPhotonMcPartObj *mymc = static_cast<AliPhotonMcPartObj*>(fMyMcParts->At(imy));
33d8da67 647 if(!mymc)
648 continue;
649 if(mymc->fLabel == iTrack)
650 checkIfAlreadySaved = true;
651 }
652 if(!checkIfAlreadySaved)
653 FillMcPart(mcP, ipart++, iTrack);
654 for(Int_t id=mcP->GetFirstDaughter(); id <= mcP->GetLastDaughter(); id++){
655 if(id<=mcP->GetMother(0))
656 continue;
657 if(id<0 || id>nTracks)
658 continue;
659 TParticle *mcD = static_cast<TParticle*>(fStack->Particle(id));
660 if(!mcD)
661 continue;
662 FillMcPart(mcD, ipart++, id);
663 }
664 }
665// cout<<"++++Get - MC done!\n";
666
667}
668//________________________________________________________________________
669void AliAnalysisTaskEMCALPhoton::FillMcPart(TParticle *mcP, Int_t ipart, Int_t itrack)
670{
671 if(!mcP)
672 return;
673 TVector3 vmcv(mcP->Vx(),mcP->Vy(), mcP->Vz());
b2d49404 674 AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(ipart));
33d8da67 675 mcp->fLabel = itrack ;
676 mcp->fPdg = mcP->GetPdgCode() ;
677 mcp->fPt = mcP->Pt() ;
678 mcp->fEta = mcP->Eta() ;
679 mcp->fPhi = mcP->Phi() ;
680 mcp->fVR = vmcv.Perp();
681 if(vmcv.Perp()>0){
682 mcp->fVEta = vmcv.Eta() ;
683 mcp->fVPhi = vmcv.Phi() ;
684 }
685 mcp->fMother = mcP->GetMother(0) ;
686}
687//________________________________________________________________________
688Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
689{
690 // Compute isolation based on tracks.
691
692 Double_t trkIsolation = 0;
693 Double_t rad2 = radius*radius;
694 Int_t ntrks = fSelPrimTracks->GetEntries();
695 for(Int_t j = 0; j<ntrks; ++j) {
696 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
697 if (!track)
698 continue;
699 if (track->Pt()<pt)
700 continue;
701
702 Float_t eta = track->Eta();
703 Float_t phi = track->Phi();
704 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
705 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
706 if(dist>rad2)
707 continue;
708 trkIsolation += track->Pt();
709 }
710 return trkIsolation;
711}
712//________________________________________________________________________
713Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const
714{
715 if(!fSelPrimTracks)
716 return 0;
717 Int_t ntracks = fSelPrimTracks->GetEntries();
718 Double_t loweta = eta - R;
719 Double_t upeta = eta + R;
720 Double_t upphi = phi + R;
721 Double_t lowphi = phi - R;
722 Double_t et = 0;
723 for(Int_t itr=0; itr<ntracks; itr++){
724 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(itr));
725 if(!track)
726 continue;
727 if(track->Pt()<minpt)
728 continue;
729 if((track->Eta() < upeta) && (track->Eta() > loweta))
730 continue;
731 if((track->Phi() > upphi) || (track->Phi() < lowphi))
732 continue;
733 et+=track->Pt();
734 }
735 return et;
736}
737//________________________________________________________________________
738Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
739{
740 // Calculate the energy of cross cells around the leading cell.
741
742 AliVCaloCells *cells = 0;
743 cells = fESD->GetEMCALCells();
744 if (!cells)
745 return 0;
746
747
748 if (!fGeom)
749 return 0;
750
751 Int_t iSupMod = -1;
752 Int_t iTower = -1;
753 Int_t iIphi = -1;
754 Int_t iIeta = -1;
755 Int_t iphi = -1;
756 Int_t ieta = -1;
757 Int_t iphis = -1;
758 Int_t ietas = -1;
759
760 Double_t crossEnergy = 0;
761
762 fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
763 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
764
765 Int_t ncells = cluster->GetNCells();
766 for (Int_t i=0; i<ncells; i++) {
767 Int_t cellAbsId = cluster->GetCellAbsId(i);
768 fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
769 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
770 Int_t aphidiff = TMath::Abs(iphi-iphis);
771 if (aphidiff>1)
772 continue;
773 Int_t aetadiff = TMath::Abs(ieta-ietas);
774 if (aetadiff>1)
775 continue;
776 if ( (aphidiff==1 && aetadiff==0) ||
777 (aphidiff==0 && aetadiff==1) ) {
778 crossEnergy += cells->GetCellAmplitude(cellAbsId);
779 }
780 }
781
782 return crossEnergy;
783}
784
785
786
787//________________________________________________________________________
788Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
789{
790 // Get maximum energy of attached cell.
791
792 id = -1;
793
794 AliVCaloCells *cells = 0;
795 cells = fESD->GetEMCALCells();
796 if (!cells)
797 return 0;
798
799 Double_t maxe = 0;
800 Int_t ncells = cluster->GetNCells();
801 for (Int_t i=0; i<ncells; i++) {
802 Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
803 if (e>maxe) {
804 maxe = e;
805 id = cluster->GetCellAbsId(i);
806 }
807 }
808 return maxe;
809}
810//________________________________________________________________________
811void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *)
812{
813 // Draw result to the screen
814 // Called once at the end of the query
815 if(fIsGrid)
816 return;
817 if (fTree) {
818 TFile *f = OpenFile(1);
819 TDirectory::TContext context(f);
820 if (f)
821 fTree->Write();
822 }
823
824
825}