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