]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALTasks/AliAnalysisTaskEMCALPhoton.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
changes in proton-lambda HBT train config macros
[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") )
772f286b 266 isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny);//kEMCEGA);
33d8da67 267
268 }
c5ac19f6 269 if(fIsMC){
270 isSelected = kTRUE;
271 }
4a4936e3 272 if(!isSelected)
273 return;
33d8da67 274
275 // Post output data.
276 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
277 if (!fESD) {
278 printf("ERROR: fESD not available, returning...\n");
279 return;
280 }
281
33d8da67 282 AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex();
283 if(!pv)
284 return;
285 if(TMath::Abs(pv->GetZ())>15)
286 return;
33d8da67 287
c5ac19f6 288 TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree();
289 TFile *inpfile = (TFile*)tree->GetCurrentFile();
290
33d8da67 291 // Track loop to fill a pT spectrum
292 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
293 AliESDtrack* track = fESD->GetTrack(iTracks);
294 if (!track)
295 continue;
296
297
298 if (fTrCuts && fTrCuts->IsSelected(track)){
299 fSelTracks->Add(track);
300 fDedxPAll->Fill(track->P(), track->GetTPCsignal());
301 }
302 if (fPrTrCuts && fPrTrCuts->IsSelected(track))
303 fSelPrimTracks->Add(track);
33d8da67 304 } //track loop
f1c66719 305
c5ac19f6 306 fHeader->fInputFileName = inpfile->GetName();
33d8da67 307 fHeader->fTrClassMask = fESD->GetHeader()->GetTriggerMask();
308 fHeader->fTrCluster = fESD->GetHeader()->GetTriggerCluster();
309 AliCentrality *cent = InputEvent()->GetCentrality();
310 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
311 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
312 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
313 if(!fIsTrain){
314 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
315 if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3)
316 break;
317 fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod);
318 }
319 }
320 fESD->GetEMCALClusters(fCaloClusters);
321 fHeader->fNClus = fCaloClusters->GetEntries();
322 fEMCalCells = fESD->GetEMCALCells();
323 fHeader->fNCells = fEMCalCells->GetNumberOfCells();
b8b781cb 324 AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts();
325 fHeader->fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5);
326
33d8da67 327 fMCEvent = MCEvent();
80c98845 328 if(fMCEvent){
33d8da67 329 fStack = (AliStack*)fMCEvent->Stack();
80c98845 330 fHeader->fNMcParts = fStack->GetNtrack();
331 }
33d8da67 332
333
334 FindConversions();
33d8da67 335 FillMyCells();
33d8da67 336 FillMyClusters();
337 if(fCaloClustersNew)
338 FillMyAltClusters();
3b37c011 339 FillIsoTracks();
33d8da67 340 if(fIsMC)
341 GetMcParts();
342
343 fTree->Fill();
344 fSelTracks->Clear();
345 fSelPrimTracks->Clear();
346 fPhotConvArray->Clear();
347 fMyClusts->Clear();
348 if(fMyAltClusts)
349 fMyAltClusts->Clear();
350 fMyCells->Clear();
351 fMyTracks->Clear();
352 if(fMyMcParts)
353 fMyMcParts->Clear();
354 fCaloClusters->Clear();
355 if(fCaloClustersNew)
356 fCaloClustersNew->Clear();
357 PostData(1, fOutputList);
735294e8 358 PostData(2, fTree);
33d8da67 359}
360
361//________________________________________________________________________
362void AliAnalysisTaskEMCALPhoton::FindConversions()
363{
f1c66719 364 // Find conversion.
365
33d8da67 366 if(!fTrCuts)
367 return;
368 Int_t iconv = 0;
369 Int_t nV0Orig = fESD->GetNumberOfV0s();
4a4936e3 370 Int_t nV0s = nV0Orig;
33d8da67 371 Int_t ntracks = fESD->GetNumberOfTracks();
4a4936e3 372 if(!fPeriod.Contains("11h") && !fPeriod.Contains("10h")){
373 fESD->ResetV0s();
374 AliV0vertexer lV0vtxer;
375 lV0vtxer.Tracks2V0vertices(fESD);
376 nV0s = fESD->GetNumberOfV0s();
377 }
33d8da67 378 fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s);
379 for(Int_t iv0=0; iv0<nV0s; iv0++){
380 AliESDv0 * v0 = fESD->GetV0(iv0);
381 if(!v0)
382 continue;
383 Double_t vpos[3];
384 fInvMassV0->Fill(v0->GetEffMass());
385 v0->GetXYZ(vpos[0], vpos[1], vpos[2]);
386 Int_t ipos = v0->GetPindex();
387 Int_t ineg = v0->GetNindex();
388 if(ipos<0 || ipos> ntracks)
389 continue;
390 if(ineg<0 || ineg> ntracks)
391 continue;
392 AliESDtrack *pos = static_cast<AliESDtrack*>(fESD->GetTrack(ipos));
393 AliESDtrack *neg = static_cast<AliESDtrack*>(fESD->GetTrack(ineg));
394 if(!pos || !neg)
395 continue;
396 if (!fTrCuts->IsSelected(pos) || !fTrCuts->IsSelected(neg))
397 continue;
398 /*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65)
399 continue;*/
400 const AliExternalTrackParam * paramPos = v0->GetParamP() ;
401 const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
402 if(!paramPos || !paramNeg)
403 continue;
404 if(pos->GetSign() <0){//change tracks
405 pos=neg ;
406 neg=fESD->GetTrack(v0->GetPindex()) ;
407 paramPos=paramNeg ;
408 paramNeg=v0->GetParamP() ;
409 }
410 AliKFParticle negKF(*paramNeg,11);
411 AliKFParticle posKF(*paramPos,-11);
412 AliKFParticle photKF(negKF,posKF) ;
33d8da67 413
414 if( neg->GetKinkIndex(0) > 0 || pos->GetKinkIndex(0) > 0)
415 continue ;
33d8da67 416
417 if(pos->GetSign() == neg->GetSign()){
418 fInvMassV0SS->Fill(photKF.GetMass());
419 continue;
420 }
421 fInvMassV0KF->Fill(photKF.GetMass());
422 TLorentzVector photLV(photKF.GetPx(), photKF.GetPy(),photKF.GetPz(), photKF.GetE());
423 if(photLV.M()>0.05 || photLV.M()<0)
424 continue;
425 fConversionVtxXY->Fill(vpos[0], vpos[1]);//negKF.GetX(), negKF.GetY());
426 Double_t convPhi = TMath::ATan2(photLV.Py(),photLV.Px());
427 if(convPhi<0)
428 convPhi+=TMath::TwoPi();
429 TVector3 vecpos(vpos);
430 Double_t v0Phi = 0;
431 if(vecpos.Perp()>0)
432 vecpos.Phi();
433 if(v0Phi<0)
434 v0Phi+=TMath::TwoPi();
b2d49404 435 AliPhotonConvObj *myconv = static_cast<AliPhotonConvObj*>(fPhotConvArray->New(iconv++));
33d8da67 436 myconv->fPt = photLV.Pt();
437 myconv->fEta = photLV.Eta();
438 myconv->fPhi = convPhi;
439 myconv->fVR = vecpos.Perp();
440 if(vecpos.Perp()>0)
441 myconv->fVEta = vecpos.Eta();
442 myconv->fVPhi = v0Phi;
443 myconv->fMass = photLV.M();
444 myconv->fMcLabel = -3; //WARNING: include the correct labeling
445 //negative daughter
446 myconv->fNegPt = negKF.GetPt();
447 myconv->fNegEta = negKF.GetEta();
448 Double_t trackPhi = negKF.GetPhi();
449 if(trackPhi<0)
450 trackPhi+=TMath::TwoPi();
451 myconv->fNegPhi = trackPhi;
452 myconv->fNegDedx = neg->GetTPCsignal();
453 myconv->fNegMcLabel = neg->GetLabel();
454 //negative daughter
455 myconv->fPosPt = posKF.GetPt();
456 myconv->fPosEta = posKF.GetEta();
457 trackPhi = posKF.GetPhi();
458 if(trackPhi<0)
459 trackPhi+=TMath::TwoPi();
460 myconv->fPosPhi = trackPhi;
461 myconv->fPosDedx = pos->GetTPCsignal();
462 myconv->fPosMcLabel = pos->GetLabel();
463 }
464 return;
465}
f1c66719 466
33d8da67 467//________________________________________________________________________
468void AliAnalysisTaskEMCALPhoton::FillMyCells()
469{
f1c66719 470 // Fill cells.
471
472 if (!fEMCalCells)
33d8da67 473 return;
33d8da67 474 Int_t ncells = fEMCalCells->GetNumberOfCells();
475 Int_t mcel = 0;
476 for(Int_t icell = 0; icell<ncells; icell++){
477 Int_t absID = TMath::Abs(fEMCalCells->GetCellNumber(icell));
b2d49404 478 AliPhotonCellObj *mycell = static_cast<AliPhotonCellObj*>(fMyCells->New(mcel++));
33d8da67 479 Float_t eta=-1, phi=-1;
480 if(!fGeom){
40b40be6 481 std::cout<<"+++fGeom not found! MyCells branch will not be filled for this event!+++\n";
33d8da67 482 return;
483 }
484 if(!fGeom)
485 return;
40b40be6 486 /*if(!fIsMC)*/fGeom->EtaPhiFromIndex(absID,eta,phi);
33d8da67 487 Float_t theta = 2*TMath::ATan(TMath::Exp(-eta));
488 mycell->fAbsID = absID;
489 mycell->fE = fEMCalCells->GetCellAmplitude(absID);
490 mycell->fEt = fEMCalCells->GetCellAmplitude(absID)*TMath::Sin(theta);
491 mycell->fEta = eta;
492 mycell->fPhi = phi;
493 mycell->fTime = fEMCalCells->GetCellTime(absID);
494 }
33d8da67 495}
f1c66719 496
33d8da67 497//________________________________________________________________________
498void AliAnalysisTaskEMCALPhoton::FillMyClusters()
499{
f1c66719 500 // Fill clusters.
501
502 if (!fCaloClusters)
33d8da67 503 return;
504 Int_t nclus = fCaloClusters->GetEntries();
505 Int_t mcl = 0;
506 for(Int_t ic=0; ic < nclus; ic++){
507 AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
508 if(!clus)
509 continue;
510 if(!clus->IsEMCAL())
511 continue;
512 if(clus->E() < fClusThresh)
513 continue;
514 Float_t pos[3];
515 clus->GetPosition(pos);
516 TVector3 cpos(pos);
517 TString cellsAbsID;
b2d49404 518 AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++));
33d8da67 519 myclus->fE = clus->E();
520 myclus->fEt = clus->E()*TMath::Sin(cpos.Theta());
521 myclus->fR = cpos.Perp();
522 myclus->fEta = cpos.Eta();
523 myclus->fPhi = cpos.Phi();
524 if(cpos.Phi()<0){
525 myclus->fPhi+=TMath::TwoPi();
526 }
527 myclus->fN = clus->GetNCells();
528 Short_t id = -1;
529 myclus->fEmax = GetMaxCellEnergy( clus, id);
530 myclus->fIdmax = id;
531 myclus->fTmax = fEMCalCells->GetCellTime(id);
532 myclus->fEcross = GetCrossEnergy( clus, id);
533 myclus->fDisp = clus->GetDispersion();
534 myclus->fM20 = clus->GetM20();
535 myclus->fM02 = clus->GetM02();
536 myclus->fTrDEta = clus->GetTrackDz();
537 myclus->fTrDPhi = clus->GetTrackDx();
538 myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
539 myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
540 myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
541 myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
542 myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.);
543 myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.);
544 myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.);
545 myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.);
546 for(Int_t icell=0;icell<myclus->fN;icell++){
547 Int_t absID = clus->GetCellAbsId(icell);
548 cellsAbsID.Append(Form("%d",absID));
549 cellsAbsID.Append(";");
550 }
551 myclus->fCellsAbsId = cellsAbsID;
552 myclus->fMcLabel = clus->GetLabel();
553 Int_t matchIndex = clus->GetTrackMatchedIndex();
554 if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
555 myclus->fTrEp = -1;
556 continue;
557 }
558 AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
559 if(!track){
560 myclus->fTrEp = -1;
561 continue;
562 }
563 if(!fPrTrCuts){
564 myclus->fTrEp = -1;
565 continue;
566 }
567 if(!fPrTrCuts->IsSelected(track)){
568 myclus->fTrEp = -1;
569 continue;
570 }
571 myclus->fTrEp = clus->E()/track->P();
572 }
573
574}
575//________________________________________________________________________
576void AliAnalysisTaskEMCALPhoton::FillMyAltClusters()
577{
f1c66719 578 // Fill clusters.
579
33d8da67 580 if(!fCaloClustersNew)
581 return;
582 Int_t nclus = fCaloClustersNew->GetEntries();
583 Int_t mcl = 0;
584 for(Int_t ic=0; ic < nclus; ic++){
585 AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClustersNew->At(ic));
586 if(!clus)
587 continue;
588 if(!clus->IsEMCAL())
589 continue;
590 if(clus->E() < fClusThresh)
591 continue;
592 Float_t pos[3];
593 clus->GetPosition(pos);
594 TVector3 cpos(pos);
595 TString cellsAbsID;
b2d49404 596 AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++));
33d8da67 597 myclus->fE = clus->E();
598 myclus->fEt = clus->E()*TMath::Sin(cpos.Theta());
599 myclus->fR = cpos.Perp();
600 myclus->fEta = cpos.Eta();
601 myclus->fPhi = cpos.Phi();
602 if(cpos.Phi()<0){
603 myclus->fPhi+=TMath::TwoPi();
604 }
605 myclus->fN = clus->GetNCells();
606 myclus->fDisp = clus->GetDispersion();
607 myclus->fM20 = clus->GetM20();
608 myclus->fM02 = clus->GetM02();
609 myclus->fTrDEta = clus->GetTrackDz();
610 myclus->fTrDPhi = clus->GetTrackDx();
611 myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
612 myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
613 myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
614 myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
615 for(Int_t icell=0;icell<myclus->fN;icell++){
616 Int_t absID = clus->GetCellAbsId(icell);
617 cellsAbsID.Append(Form("%d",absID));
618 cellsAbsID.Append(";");
619 }
620 myclus->fCellsAbsId = cellsAbsID;
621 myclus->fMcLabel = clus->GetLabel();
622 Int_t matchIndex = clus->GetTrackMatchedIndex();
623 if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
624 myclus->fTrEp = -1;
625 continue;
626 }
627 AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
628 if(!track){
629 myclus->fTrEp = -1;
630 continue;
631 }
632 if(!fPrTrCuts){
633 myclus->fTrEp = -1;
634 continue;
635 }
636 if(!fPrTrCuts->IsSelected(track)){
637 myclus->fTrEp = -1;
638 continue;
639 }
640 myclus->fTrEp = clus->E()/track->P();
641 }
642
643}
644//________________________________________________________________________
3b37c011 645void AliAnalysisTaskEMCALPhoton::FillIsoTracks()
33d8da67 646{
f1c66719 647 // Fill high pt tracks.
648
33d8da67 649 if(!fSelPrimTracks)
650 return;
651 Int_t ntracks = fSelPrimTracks->GetEntries();
652 Int_t imtr = 0;
653 for(Int_t it=0;it<ntracks; it++){
654 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(it));
655 if(!track)
656 continue;
3b37c011 657 if(track->Phi()<1.0 || track->Phi()>3.55)
658 continue;
659 if(TMath::Abs(track->Eta())>1.1)
660 continue;
a72860a2 661 /*if(track->Pt()<3)
662 continue;*/
663 AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
33d8da67 664 mtr->fPt = track->Pt();
665 mtr->fEta = track->Eta();
666 mtr->fPhi = track->Phi();
3b37c011 667 mtr->fCharge = track->Charge();
33d8da67 668 mtr->fDedx = track->GetTPCsignal();
669 mtr->fMcLabel = track->GetLabel();
670 }
671}
f1c66719 672
33d8da67 673//________________________________________________________________________
674void AliAnalysisTaskEMCALPhoton::GetMcParts()
675{
f1c66719 676 // Get MC particles.
677
33d8da67 678 if (!fStack)
679 return;
33d8da67 680
681 const AliVVertex *evtVtx = fMCEvent->GetPrimaryVertex();
682 if (!evtVtx)
683 return;
33d8da67 684 Int_t nTracks = fStack->GetNtrack();
80c98845 685 AliPhotonMcPartObj *mcp = 0;
33d8da67 686 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
687 TParticle *mcP = static_cast<TParticle*>(fStack->Particle(iTrack));
80c98845 688 if (!mcP){
689 mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
df3bb571 690 mcp->fLabel = iTrack;
33d8da67 691 continue;
80c98845 692 }
33d8da67 693 // primary particle
694 Double_t dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
695 (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
696 (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
80c98845 697 if(dR > 0.5){
698 mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
df3bb571 699 mcp->fLabel = iTrack;
33d8da67 700 continue;
80c98845 701 }
33d8da67 702
703 // kinematic cuts
704 Double_t pt = mcP->Pt() ;
80c98845 705 if (pt<0.5){
706 mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
df3bb571 707 mcp->fLabel = iTrack;
33d8da67 708 continue;
80c98845 709 }
33d8da67 710 Double_t eta = mcP->Eta();
80c98845 711 if (TMath::Abs(eta)>0.7){
712 mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
df3bb571 713 mcp->fLabel = iTrack;
33d8da67 714 continue;
80c98845 715 }
33d8da67 716 Double_t phi = mcP->Phi();
80c98845 717 if (phi<1.0||phi>3.3){
718 mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
e68b22e2 719 mcp->fLabel = iTrack;
33d8da67 720 continue;
80c98845 721 }
33d8da67 722 // pion or eta meson or direct photon
723 if(mcP->GetPdgCode() == 111) {
724 } else if(mcP->GetPdgCode() == 221) {
725 } else if(mcP->GetPdgCode() == 22 ) {
80c98845 726 } else {
727 mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(iTrack));
df3bb571 728 mcp->fLabel = iTrack;
33d8da67 729 continue;
80c98845 730 }
731
33d8da67 732 bool checkIfAlreadySaved = false;
733 for(Int_t imy=0;imy<fMyMcParts->GetEntries();imy++){
b2d49404 734 AliPhotonMcPartObj *mymc = static_cast<AliPhotonMcPartObj*>(fMyMcParts->At(imy));
33d8da67 735 if(!mymc)
736 continue;
80c98845 737 if(imy == iTrack)
33d8da67 738 checkIfAlreadySaved = true;
739 }
80c98845 740 if(!checkIfAlreadySaved){
741 FillMcPart(mcP, iTrack);
742 }
33d8da67 743 for(Int_t id=mcP->GetFirstDaughter(); id <= mcP->GetLastDaughter(); id++){
744 if(id<=mcP->GetMother(0))
745 continue;
746 if(id<0 || id>nTracks)
747 continue;
748 TParticle *mcD = static_cast<TParticle*>(fStack->Particle(id));
749 if(!mcD)
750 continue;
80c98845 751 FillMcPart(mcD, id);
752 }
33d8da67 753 }
33d8da67 754}
f1c66719 755
33d8da67 756//________________________________________________________________________
80c98845 757void AliAnalysisTaskEMCALPhoton::FillMcPart(TParticle *mcP, Int_t itrack)
33d8da67 758{
f1c66719 759 // Fill MC particles.
760
33d8da67 761 if(!mcP)
762 return;
763 TVector3 vmcv(mcP->Vx(),mcP->Vy(), mcP->Vz());
80c98845 764 AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(itrack));
33d8da67 765 mcp->fLabel = itrack ;
766 mcp->fPdg = mcP->GetPdgCode() ;
767 mcp->fPt = mcP->Pt() ;
768 mcp->fEta = mcP->Eta() ;
769 mcp->fPhi = mcP->Phi() ;
770 mcp->fVR = vmcv.Perp();
771 if(vmcv.Perp()>0){
772 mcp->fVEta = vmcv.Eta() ;
773 mcp->fVPhi = vmcv.Phi() ;
774 }
775 mcp->fMother = mcP->GetMother(0) ;
80c98845 776 mcp->fFirstD = mcP->GetFirstDaughter() ;
777 mcp->fLastD = mcP->GetLastDaughter() ;
778 mcp->fStatus = mcP->GetStatusCode();
779 mcp->fIso = AliAnalysisTaskEMCALPhoton::GetMcIsolation(mcP, itrack, 0.4 , 0.2);
780}
781//________________________________________________________________________
782Double_t AliAnalysisTaskEMCALPhoton::GetMcIsolation(TParticle *mcP, Int_t itrack, Double_t radius, Double_t pt) const
783{
784 if (!fStack)
785 return -1;
786 if(!mcP)
787 return -1;
788 if(itrack<6 || itrack>8)
789 return -1;
790 if(mcP->GetPdgCode()!=22)
791 return -1;
792 Double_t sumpt=0;
793 Float_t eta = mcP->Eta();
794 Float_t phi = mcP->Phi();
795 Float_t dR;
796 Int_t nparts = fStack->GetNtrack();
797 for(Int_t ip = 0; ip<nparts; ip++){
798 TParticle *mcisop = static_cast<TParticle*>(fStack->Particle(ip));
799 if(!mcisop)
800 continue;
801 if(ip==itrack)
802 continue;
803 if(mcisop->GetStatusCode()!=1)
804 continue;
805 if(mcisop->Pt()<pt)
806 continue;
807 TVector3 vmcv(mcisop->Vx(),mcisop->Vy(), mcisop->Vz());
808 if(vmcv.Perp()>1)
809 continue;
810 dR = TMath::Sqrt((phi-mcisop->Phi())*(phi-mcisop->Phi())+(eta-mcisop->Eta())*(eta-mcisop->Eta()));
811 if(dR>radius)
812 continue;
813 sumpt += mcisop->Pt();
814 }
815 return sumpt;
33d8da67 816}
f1c66719 817
33d8da67 818//________________________________________________________________________
819Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
820{
821 // Compute isolation based on tracks.
822
823 Double_t trkIsolation = 0;
824 Double_t rad2 = radius*radius;
825 Int_t ntrks = fSelPrimTracks->GetEntries();
826 for(Int_t j = 0; j<ntrks; ++j) {
827 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
828 if (!track)
829 continue;
830 if (track->Pt()<pt)
831 continue;
832
833 Float_t eta = track->Eta();
834 Float_t phi = track->Phi();
835 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
836 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
837 if(dist>rad2)
838 continue;
839 trkIsolation += track->Pt();
840 }
841 return trkIsolation;
842}
f1c66719 843
33d8da67 844//________________________________________________________________________
845Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const
846{
f1c66719 847 // Get phi band.
848
33d8da67 849 if(!fSelPrimTracks)
850 return 0;
851 Int_t ntracks = fSelPrimTracks->GetEntries();
852 Double_t loweta = eta - R;
853 Double_t upeta = eta + R;
854 Double_t upphi = phi + R;
855 Double_t lowphi = phi - R;
856 Double_t et = 0;
857 for(Int_t itr=0; itr<ntracks; itr++){
858 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(itr));
859 if(!track)
860 continue;
861 if(track->Pt()<minpt)
862 continue;
863 if((track->Eta() < upeta) && (track->Eta() > loweta))
864 continue;
865 if((track->Phi() > upphi) || (track->Phi() < lowphi))
866 continue;
867 et+=track->Pt();
868 }
869 return et;
870}
f1c66719 871
33d8da67 872//________________________________________________________________________
873Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
874{
875 // Calculate the energy of cross cells around the leading cell.
876
877 AliVCaloCells *cells = 0;
878 cells = fESD->GetEMCALCells();
879 if (!cells)
880 return 0;
881
33d8da67 882 if (!fGeom)
883 return 0;
884
885 Int_t iSupMod = -1;
886 Int_t iTower = -1;
887 Int_t iIphi = -1;
888 Int_t iIeta = -1;
889 Int_t iphi = -1;
890 Int_t ieta = -1;
891 Int_t iphis = -1;
892 Int_t ietas = -1;
893
894 Double_t crossEnergy = 0;
895
896 fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
897 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
898
899 Int_t ncells = cluster->GetNCells();
900 for (Int_t i=0; i<ncells; i++) {
901 Int_t cellAbsId = cluster->GetCellAbsId(i);
902 fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
903 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
904 Int_t aphidiff = TMath::Abs(iphi-iphis);
905 if (aphidiff>1)
906 continue;
907 Int_t aetadiff = TMath::Abs(ieta-ietas);
908 if (aetadiff>1)
909 continue;
910 if ( (aphidiff==1 && aetadiff==0) ||
911 (aphidiff==0 && aetadiff==1) ) {
912 crossEnergy += cells->GetCellAmplitude(cellAbsId);
913 }
914 }
915
916 return crossEnergy;
917}
918
33d8da67 919//________________________________________________________________________
920Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
921{
922 // Get maximum energy of attached cell.
923
924 id = -1;
925
926 AliVCaloCells *cells = 0;
927 cells = fESD->GetEMCALCells();
928 if (!cells)
929 return 0;
930
931 Double_t maxe = 0;
932 Int_t ncells = cluster->GetNCells();
933 for (Int_t i=0; i<ncells; i++) {
934 Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
935 if (e>maxe) {
936 maxe = e;
937 id = cluster->GetCellAbsId(i);
938 }
939 }
940 return maxe;
941}
f1c66719 942
33d8da67 943//________________________________________________________________________
944void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *)
945{
33d8da67 946 // Called once at the end of the query
ec681cf3 947/* if(fIsGrid)
948 return;*/
f1c66719 949 if (fTree) {
c305511b 950 printf("***tree %s being saved***\n",fTree->GetName());
55dcd0ae 951 TFile *f = OpenFile(2);
33d8da67 952 TDirectory::TContext context(f);
953 if (f)
954 fTree->Write();
955 }
33d8da67 956}