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