disabling redo V0 vertexer in PbPb and warning fix
[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 }
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();
328 if(fMCEvent)
329 fStack = (AliStack*)fMCEvent->Stack();
330
331
332 FindConversions();
33d8da67 333 FillMyCells();
33d8da67 334 FillMyClusters();
335 if(fCaloClustersNew)
336 FillMyAltClusters();
3b37c011 337 FillIsoTracks();
33d8da67 338 if(fIsMC)
339 GetMcParts();
340
341 fTree->Fill();
342 fSelTracks->Clear();
343 fSelPrimTracks->Clear();
344 fPhotConvArray->Clear();
345 fMyClusts->Clear();
346 if(fMyAltClusts)
347 fMyAltClusts->Clear();
348 fMyCells->Clear();
349 fMyTracks->Clear();
350 if(fMyMcParts)
351 fMyMcParts->Clear();
352 fCaloClusters->Clear();
353 if(fCaloClustersNew)
354 fCaloClustersNew->Clear();
355 PostData(1, fOutputList);
735294e8 356 PostData(2, fTree);
33d8da67 357}
358
359//________________________________________________________________________
360void AliAnalysisTaskEMCALPhoton::FindConversions()
361{
f1c66719 362 // Find conversion.
363
33d8da67 364 if(!fTrCuts)
365 return;
366 Int_t iconv = 0;
367 Int_t nV0Orig = fESD->GetNumberOfV0s();
4a4936e3 368 Int_t nV0s = nV0Orig;
33d8da67 369 Int_t ntracks = fESD->GetNumberOfTracks();
4a4936e3 370 if(!fPeriod.Contains("11h") && !fPeriod.Contains("10h")){
371 fESD->ResetV0s();
372 AliV0vertexer lV0vtxer;
373 lV0vtxer.Tracks2V0vertices(fESD);
374 nV0s = fESD->GetNumberOfV0s();
375 }
33d8da67 376 fNV0sBefAndAftRerun->Fill(nV0Orig, nV0s);
377 for(Int_t iv0=0; iv0<nV0s; iv0++){
378 AliESDv0 * v0 = fESD->GetV0(iv0);
379 if(!v0)
380 continue;
381 Double_t vpos[3];
382 fInvMassV0->Fill(v0->GetEffMass());
383 v0->GetXYZ(vpos[0], vpos[1], vpos[2]);
384 Int_t ipos = v0->GetPindex();
385 Int_t ineg = v0->GetNindex();
386 if(ipos<0 || ipos> ntracks)
387 continue;
388 if(ineg<0 || ineg> ntracks)
389 continue;
390 AliESDtrack *pos = static_cast<AliESDtrack*>(fESD->GetTrack(ipos));
391 AliESDtrack *neg = static_cast<AliESDtrack*>(fESD->GetTrack(ineg));
392 if(!pos || !neg)
393 continue;
394 if (!fTrCuts->IsSelected(pos) || !fTrCuts->IsSelected(neg))
395 continue;
396 /*if(pos->GetTPCsignal()<65 || neg->GetTPCsignal()<65)
397 continue;*/
398 const AliExternalTrackParam * paramPos = v0->GetParamP() ;
399 const AliExternalTrackParam * paramNeg = v0->GetParamN() ;
400 if(!paramPos || !paramNeg)
401 continue;
402 if(pos->GetSign() <0){//change tracks
403 pos=neg ;
404 neg=fESD->GetTrack(v0->GetPindex()) ;
405 paramPos=paramNeg ;
406 paramNeg=v0->GetParamP() ;
407 }
408 AliKFParticle negKF(*paramNeg,11);
409 AliKFParticle posKF(*paramPos,-11);
410 AliKFParticle photKF(negKF,posKF) ;
33d8da67 411
412 if( neg->GetKinkIndex(0) > 0 || pos->GetKinkIndex(0) > 0)
413 continue ;
33d8da67 414
415 if(pos->GetSign() == neg->GetSign()){
416 fInvMassV0SS->Fill(photKF.GetMass());
417 continue;
418 }
419 fInvMassV0KF->Fill(photKF.GetMass());
420 TLorentzVector photLV(photKF.GetPx(), photKF.GetPy(),photKF.GetPz(), photKF.GetE());
421 if(photLV.M()>0.05 || photLV.M()<0)
422 continue;
423 fConversionVtxXY->Fill(vpos[0], vpos[1]);//negKF.GetX(), negKF.GetY());
424 Double_t convPhi = TMath::ATan2(photLV.Py(),photLV.Px());
425 if(convPhi<0)
426 convPhi+=TMath::TwoPi();
427 TVector3 vecpos(vpos);
428 Double_t v0Phi = 0;
429 if(vecpos.Perp()>0)
430 vecpos.Phi();
431 if(v0Phi<0)
432 v0Phi+=TMath::TwoPi();
b2d49404 433 AliPhotonConvObj *myconv = static_cast<AliPhotonConvObj*>(fPhotConvArray->New(iconv++));
33d8da67 434 myconv->fPt = photLV.Pt();
435 myconv->fEta = photLV.Eta();
436 myconv->fPhi = convPhi;
437 myconv->fVR = vecpos.Perp();
438 if(vecpos.Perp()>0)
439 myconv->fVEta = vecpos.Eta();
440 myconv->fVPhi = v0Phi;
441 myconv->fMass = photLV.M();
442 myconv->fMcLabel = -3; //WARNING: include the correct labeling
443 //negative daughter
444 myconv->fNegPt = negKF.GetPt();
445 myconv->fNegEta = negKF.GetEta();
446 Double_t trackPhi = negKF.GetPhi();
447 if(trackPhi<0)
448 trackPhi+=TMath::TwoPi();
449 myconv->fNegPhi = trackPhi;
450 myconv->fNegDedx = neg->GetTPCsignal();
451 myconv->fNegMcLabel = neg->GetLabel();
452 //negative daughter
453 myconv->fPosPt = posKF.GetPt();
454 myconv->fPosEta = posKF.GetEta();
455 trackPhi = posKF.GetPhi();
456 if(trackPhi<0)
457 trackPhi+=TMath::TwoPi();
458 myconv->fPosPhi = trackPhi;
459 myconv->fPosDedx = pos->GetTPCsignal();
460 myconv->fPosMcLabel = pos->GetLabel();
461 }
462 return;
463}
f1c66719 464
33d8da67 465//________________________________________________________________________
466void AliAnalysisTaskEMCALPhoton::FillMyCells()
467{
f1c66719 468 // Fill cells.
469
470 if (!fEMCalCells)
33d8da67 471 return;
33d8da67 472 Int_t ncells = fEMCalCells->GetNumberOfCells();
473 Int_t mcel = 0;
474 for(Int_t icell = 0; icell<ncells; icell++){
475 Int_t absID = TMath::Abs(fEMCalCells->GetCellNumber(icell));
b2d49404 476 AliPhotonCellObj *mycell = static_cast<AliPhotonCellObj*>(fMyCells->New(mcel++));
33d8da67 477 Float_t eta=-1, phi=-1;
478 if(!fGeom){
40b40be6 479 std::cout<<"+++fGeom not found! MyCells branch will not be filled for this event!+++\n";
33d8da67 480 return;
481 }
482 if(!fGeom)
483 return;
40b40be6 484 /*if(!fIsMC)*/fGeom->EtaPhiFromIndex(absID,eta,phi);
33d8da67 485 Float_t theta = 2*TMath::ATan(TMath::Exp(-eta));
486 mycell->fAbsID = absID;
487 mycell->fE = fEMCalCells->GetCellAmplitude(absID);
488 mycell->fEt = fEMCalCells->GetCellAmplitude(absID)*TMath::Sin(theta);
489 mycell->fEta = eta;
490 mycell->fPhi = phi;
491 mycell->fTime = fEMCalCells->GetCellTime(absID);
492 }
33d8da67 493}
f1c66719 494
33d8da67 495//________________________________________________________________________
496void AliAnalysisTaskEMCALPhoton::FillMyClusters()
497{
f1c66719 498 // Fill clusters.
499
500 if (!fCaloClusters)
33d8da67 501 return;
502 Int_t nclus = fCaloClusters->GetEntries();
503 Int_t mcl = 0;
504 for(Int_t ic=0; ic < nclus; ic++){
505 AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClusters->At(ic));
506 if(!clus)
507 continue;
508 if(!clus->IsEMCAL())
509 continue;
510 if(clus->E() < fClusThresh)
511 continue;
512 Float_t pos[3];
513 clus->GetPosition(pos);
514 TVector3 cpos(pos);
515 TString cellsAbsID;
b2d49404 516 AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyClusts->New(mcl++));
33d8da67 517 myclus->fE = clus->E();
518 myclus->fEt = clus->E()*TMath::Sin(cpos.Theta());
519 myclus->fR = cpos.Perp();
520 myclus->fEta = cpos.Eta();
521 myclus->fPhi = cpos.Phi();
522 if(cpos.Phi()<0){
523 myclus->fPhi+=TMath::TwoPi();
524 }
525 myclus->fN = clus->GetNCells();
526 Short_t id = -1;
527 myclus->fEmax = GetMaxCellEnergy( clus, id);
528 myclus->fIdmax = id;
529 myclus->fTmax = fEMCalCells->GetCellTime(id);
530 myclus->fEcross = GetCrossEnergy( clus, id);
531 myclus->fDisp = clus->GetDispersion();
532 myclus->fM20 = clus->GetM20();
533 myclus->fM02 = clus->GetM02();
534 myclus->fTrDEta = clus->GetTrackDz();
535 myclus->fTrDPhi = clus->GetTrackDx();
536 myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
537 myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
538 myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
539 myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
540 myclus->fTrPhiBand01 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.1, 0.);
541 myclus->fTrPhiBand02 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.2, 0.);
542 myclus->fTrPhiBand03 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.3, 0.);
543 myclus->fTrPhiBand04 = GetPhiBandEt( cpos.Eta(), cpos.Phi(), 0.4, 0.);
544 for(Int_t icell=0;icell<myclus->fN;icell++){
545 Int_t absID = clus->GetCellAbsId(icell);
546 cellsAbsID.Append(Form("%d",absID));
547 cellsAbsID.Append(";");
548 }
549 myclus->fCellsAbsId = cellsAbsID;
550 myclus->fMcLabel = clus->GetLabel();
551 Int_t matchIndex = clus->GetTrackMatchedIndex();
552 if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
553 myclus->fTrEp = -1;
554 continue;
555 }
556 AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
557 if(!track){
558 myclus->fTrEp = -1;
559 continue;
560 }
561 if(!fPrTrCuts){
562 myclus->fTrEp = -1;
563 continue;
564 }
565 if(!fPrTrCuts->IsSelected(track)){
566 myclus->fTrEp = -1;
567 continue;
568 }
569 myclus->fTrEp = clus->E()/track->P();
570 }
571
572}
573//________________________________________________________________________
574void AliAnalysisTaskEMCALPhoton::FillMyAltClusters()
575{
f1c66719 576 // Fill clusters.
577
33d8da67 578 if(!fCaloClustersNew)
579 return;
580 Int_t nclus = fCaloClustersNew->GetEntries();
581 Int_t mcl = 0;
582 for(Int_t ic=0; ic < nclus; ic++){
583 AliESDCaloCluster *clus = static_cast<AliESDCaloCluster*>(fCaloClustersNew->At(ic));
584 if(!clus)
585 continue;
586 if(!clus->IsEMCAL())
587 continue;
588 if(clus->E() < fClusThresh)
589 continue;
590 Float_t pos[3];
591 clus->GetPosition(pos);
592 TVector3 cpos(pos);
593 TString cellsAbsID;
b2d49404 594 AliPhotonClusterObj *myclus = static_cast<AliPhotonClusterObj*>(fMyAltClusts->New(mcl++));
33d8da67 595 myclus->fE = clus->E();
596 myclus->fEt = clus->E()*TMath::Sin(cpos.Theta());
597 myclus->fR = cpos.Perp();
598 myclus->fEta = cpos.Eta();
599 myclus->fPhi = cpos.Phi();
600 if(cpos.Phi()<0){
601 myclus->fPhi+=TMath::TwoPi();
602 }
603 myclus->fN = clus->GetNCells();
604 myclus->fDisp = clus->GetDispersion();
605 myclus->fM20 = clus->GetM20();
606 myclus->fM02 = clus->GetM02();
607 myclus->fTrDEta = clus->GetTrackDz();
608 myclus->fTrDPhi = clus->GetTrackDx();
609 myclus->fTrIso01 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.1, 0.);
610 myclus->fTrIso02 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.2, 0.);
611 myclus->fTrIso03 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.3, 0.);
612 myclus->fTrIso04 = GetTrackIsolation( cpos.Eta(), cpos.Phi(), 0.4, 0.);
613 for(Int_t icell=0;icell<myclus->fN;icell++){
614 Int_t absID = clus->GetCellAbsId(icell);
615 cellsAbsID.Append(Form("%d",absID));
616 cellsAbsID.Append(";");
617 }
618 myclus->fCellsAbsId = cellsAbsID;
619 myclus->fMcLabel = clus->GetLabel();
620 Int_t matchIndex = clus->GetTrackMatchedIndex();
621 if(matchIndex<0 || matchIndex>fESD->GetNumberOfTracks()){
622 myclus->fTrEp = -1;
623 continue;
624 }
625 AliESDtrack* track = static_cast<AliESDtrack*>(fESD->GetTrack(matchIndex));
626 if(!track){
627 myclus->fTrEp = -1;
628 continue;
629 }
630 if(!fPrTrCuts){
631 myclus->fTrEp = -1;
632 continue;
633 }
634 if(!fPrTrCuts->IsSelected(track)){
635 myclus->fTrEp = -1;
636 continue;
637 }
638 myclus->fTrEp = clus->E()/track->P();
639 }
640
641}
642//________________________________________________________________________
3b37c011 643void AliAnalysisTaskEMCALPhoton::FillIsoTracks()
33d8da67 644{
f1c66719 645 // Fill high pt tracks.
646
33d8da67 647 if(!fSelPrimTracks)
648 return;
649 Int_t ntracks = fSelPrimTracks->GetEntries();
650 Int_t imtr = 0;
651 for(Int_t it=0;it<ntracks; it++){
652 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(it));
653 if(!track)
654 continue;
3b37c011 655 /*if(track->Pt()<3)
656 continue;*/
b2d49404 657 AliPhotonTrackObj *mtr = static_cast<AliPhotonTrackObj*>(fMyTracks->New(imtr++));
3b37c011 658 if(track->Phi()<1.0 || track->Phi()>3.55)
659 continue;
660 if(TMath::Abs(track->Eta())>1.1)
661 continue;
33d8da67 662 mtr->fPt = track->Pt();
663 mtr->fEta = track->Eta();
664 mtr->fPhi = track->Phi();
3b37c011 665 mtr->fCharge = track->Charge();
33d8da67 666 mtr->fDedx = track->GetTPCsignal();
667 mtr->fMcLabel = track->GetLabel();
668 }
669}
f1c66719 670
33d8da67 671//________________________________________________________________________
672void AliAnalysisTaskEMCALPhoton::GetMcParts()
673{
f1c66719 674 // Get MC particles.
675
33d8da67 676 if (!fStack)
677 return;
33d8da67 678
679 const AliVVertex *evtVtx = fMCEvent->GetPrimaryVertex();
680 if (!evtVtx)
681 return;
682 Int_t ipart = 0;
683 Int_t nTracks = fStack->GetNtrack();
684 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
685 TParticle *mcP = static_cast<TParticle*>(fStack->Particle(iTrack));
686 if (!mcP)
687 continue;
688 // primary particle
689 Double_t dR = TMath::Sqrt((mcP->Vx()-evtVtx->GetX())*(mcP->Vx()-evtVtx->GetX()) +
690 (mcP->Vy()-evtVtx->GetY())*(mcP->Vy()-evtVtx->GetY()) +
691 (mcP->Vz()-evtVtx->GetZ())*(mcP->Vz()-evtVtx->GetZ()));
692 if(dR > 0.5)
693 continue;
694
695 // kinematic cuts
696 Double_t pt = mcP->Pt() ;
697 if (pt<0.5)
698 continue;
699 Double_t eta = mcP->Eta();
700 if (eta<-0.7||eta>0.7)
701 continue;
702 Double_t phi = mcP->Phi();
703 if (phi<1.0||phi>3.3)
704 continue;
705 // pion or eta meson or direct photon
706 if(mcP->GetPdgCode() == 111) {
707 } else if(mcP->GetPdgCode() == 221) {
708 } else if(mcP->GetPdgCode() == 22 ) {
709 } else
710 continue;
711 bool checkIfAlreadySaved = false;
712 for(Int_t imy=0;imy<fMyMcParts->GetEntries();imy++){
b2d49404 713 AliPhotonMcPartObj *mymc = static_cast<AliPhotonMcPartObj*>(fMyMcParts->At(imy));
33d8da67 714 if(!mymc)
715 continue;
716 if(mymc->fLabel == iTrack)
717 checkIfAlreadySaved = true;
718 }
719 if(!checkIfAlreadySaved)
720 FillMcPart(mcP, ipart++, iTrack);
721 for(Int_t id=mcP->GetFirstDaughter(); id <= mcP->GetLastDaughter(); id++){
722 if(id<=mcP->GetMother(0))
723 continue;
724 if(id<0 || id>nTracks)
725 continue;
726 TParticle *mcD = static_cast<TParticle*>(fStack->Particle(id));
727 if(!mcD)
728 continue;
729 FillMcPart(mcD, ipart++, id);
730 }
731 }
33d8da67 732}
f1c66719 733
33d8da67 734//________________________________________________________________________
735void AliAnalysisTaskEMCALPhoton::FillMcPart(TParticle *mcP, Int_t ipart, Int_t itrack)
736{
f1c66719 737 // Fill MC particles.
738
33d8da67 739 if(!mcP)
740 return;
741 TVector3 vmcv(mcP->Vx(),mcP->Vy(), mcP->Vz());
b2d49404 742 AliPhotonMcPartObj *mcp = static_cast<AliPhotonMcPartObj*>(fMyMcParts->New(ipart));
33d8da67 743 mcp->fLabel = itrack ;
744 mcp->fPdg = mcP->GetPdgCode() ;
745 mcp->fPt = mcP->Pt() ;
746 mcp->fEta = mcP->Eta() ;
747 mcp->fPhi = mcP->Phi() ;
748 mcp->fVR = vmcv.Perp();
749 if(vmcv.Perp()>0){
750 mcp->fVEta = vmcv.Eta() ;
751 mcp->fVPhi = vmcv.Phi() ;
752 }
753 mcp->fMother = mcP->GetMother(0) ;
754}
f1c66719 755
33d8da67 756//________________________________________________________________________
757Double_t AliAnalysisTaskEMCALPhoton::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
758{
759 // Compute isolation based on tracks.
760
761 Double_t trkIsolation = 0;
762 Double_t rad2 = radius*radius;
763 Int_t ntrks = fSelPrimTracks->GetEntries();
764 for(Int_t j = 0; j<ntrks; ++j) {
765 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
766 if (!track)
767 continue;
768 if (track->Pt()<pt)
769 continue;
770
771 Float_t eta = track->Eta();
772 Float_t phi = track->Phi();
773 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
774 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
775 if(dist>rad2)
776 continue;
777 trkIsolation += track->Pt();
778 }
779 return trkIsolation;
780}
f1c66719 781
33d8da67 782//________________________________________________________________________
783Double_t AliAnalysisTaskEMCALPhoton::GetPhiBandEt(Double_t eta, Double_t phi, Double_t R, Double_t minpt) const
784{
f1c66719 785 // Get phi band.
786
33d8da67 787 if(!fSelPrimTracks)
788 return 0;
789 Int_t ntracks = fSelPrimTracks->GetEntries();
790 Double_t loweta = eta - R;
791 Double_t upeta = eta + R;
792 Double_t upphi = phi + R;
793 Double_t lowphi = phi - R;
794 Double_t et = 0;
795 for(Int_t itr=0; itr<ntracks; itr++){
796 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(itr));
797 if(!track)
798 continue;
799 if(track->Pt()<minpt)
800 continue;
801 if((track->Eta() < upeta) && (track->Eta() > loweta))
802 continue;
803 if((track->Phi() > upphi) || (track->Phi() < lowphi))
804 continue;
805 et+=track->Pt();
806 }
807 return et;
808}
f1c66719 809
33d8da67 810//________________________________________________________________________
811Double_t AliAnalysisTaskEMCALPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax)
812{
813 // Calculate the energy of cross cells around the leading cell.
814
815 AliVCaloCells *cells = 0;
816 cells = fESD->GetEMCALCells();
817 if (!cells)
818 return 0;
819
33d8da67 820 if (!fGeom)
821 return 0;
822
823 Int_t iSupMod = -1;
824 Int_t iTower = -1;
825 Int_t iIphi = -1;
826 Int_t iIeta = -1;
827 Int_t iphi = -1;
828 Int_t ieta = -1;
829 Int_t iphis = -1;
830 Int_t ietas = -1;
831
832 Double_t crossEnergy = 0;
833
834 fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta);
835 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas);
836
837 Int_t ncells = cluster->GetNCells();
838 for (Int_t i=0; i<ncells; i++) {
839 Int_t cellAbsId = cluster->GetCellAbsId(i);
840 fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta);
841 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta);
842 Int_t aphidiff = TMath::Abs(iphi-iphis);
843 if (aphidiff>1)
844 continue;
845 Int_t aetadiff = TMath::Abs(ieta-ietas);
846 if (aetadiff>1)
847 continue;
848 if ( (aphidiff==1 && aetadiff==0) ||
849 (aphidiff==0 && aetadiff==1) ) {
850 crossEnergy += cells->GetCellAmplitude(cellAbsId);
851 }
852 }
853
854 return crossEnergy;
855}
856
33d8da67 857//________________________________________________________________________
858Double_t AliAnalysisTaskEMCALPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
859{
860 // Get maximum energy of attached cell.
861
862 id = -1;
863
864 AliVCaloCells *cells = 0;
865 cells = fESD->GetEMCALCells();
866 if (!cells)
867 return 0;
868
869 Double_t maxe = 0;
870 Int_t ncells = cluster->GetNCells();
871 for (Int_t i=0; i<ncells; i++) {
872 Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
873 if (e>maxe) {
874 maxe = e;
875 id = cluster->GetCellAbsId(i);
876 }
877 }
878 return maxe;
879}
f1c66719 880
33d8da67 881//________________________________________________________________________
882void AliAnalysisTaskEMCALPhoton::Terminate(Option_t *)
883{
33d8da67 884 // Called once at the end of the query
ec681cf3 885/* if(fIsGrid)
886 return;*/
f1c66719 887 if (fTree) {
c305511b 888 printf("***tree %s being saved***\n",fTree->GetName());
55dcd0ae 889 TFile *f = OpenFile(2);
33d8da67 890 TDirectory::TContext context(f);
891 if (f)
892 fTree->Write();
893 }
33d8da67 894}