]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/UserTasks/EmcalTasks/AliAnalysisTaskEMCALPi0PbPb.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
sort digits
[u/mrichter/AliRoot.git] / PWG4 / UserTasks / EmcalTasks / AliAnalysisTaskEMCALPi0PbPb.cxx
CommitLineData
ea3fd2d5 1// $Id$
2
ea3fd2d5 3#include "AliAnalysisTaskEMCALPi0PbPb.h"
fa443410 4#include <TAxis.h>
717fe7de 5#include <TChain.h>
6#include <TClonesArray.h>
f5d4ab70 7#include <TFile.h>
296ea9b4 8#include <TGeoGlobalMagField.h>
717fe7de 9#include <TH1F.h>
10#include <TH2F.h>
11#include <TList.h>
12#include <TLorentzVector.h>
f5d4ab70 13#include <TNtuple.h>
296ea9b4 14#include <TProfile.h>
0ceca1f8 15#include <TRegexp.h>
b3ee6797 16#include <TString.h>
296ea9b4 17#include <TVector2.h>
717fe7de 18#include "AliAODEvent.h"
56fd6cb2 19#include "AliAODMCParticle.h"
717fe7de 20#include "AliAODVertex.h"
21#include "AliAnalysisManager.h"
22#include "AliAnalysisTaskEMCALClusterizeFast.h"
296ea9b4 23#include "AliCDBManager.h"
ea3fd2d5 24#include "AliCentrality.h"
c2fe5f0e 25#include "AliEMCALDigit.h"
2ef5608f 26#include "AliEMCALGeometry.h"
0fbe8d4f 27#include "AliEMCALRecPoint.h"
296ea9b4 28#include "AliEMCALRecoUtils.h"
3a952328 29#include "AliESDCaloTrigger.h"
717fe7de 30#include "AliESDEvent.h"
5fe1ca23 31#include "AliESDUtils.h"
717fe7de 32#include "AliESDVertex.h"
296ea9b4 33#include "AliESDtrack.h"
0ec74551 34#include "AliESDtrackCuts.h"
b6c599fe 35#include "AliEventplane.h"
296ea9b4 36#include "AliGeomManager.h"
b3ee6797 37#include "AliInputEventHandler.h"
43cfaa06 38#include "AliLog.h"
38727e64 39#include "AliMCEvent.h"
40#include "AliMCParticle.h"
296ea9b4 41#include "AliMagF.h"
5fe1ca23 42#include "AliMultiplicity.h"
38727e64 43#include "AliStack.h"
0ec74551 44#include "AliTrackerBase.h"
c0d2e1f2 45#include "AliTriggerAnalysis.h"
ea3fd2d5 46
47ClassImp(AliAnalysisTaskEMCALPi0PbPb)
717fe7de 48
4ea96211 49//________________________________________________________________________
50AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb()
51 : AliAnalysisTaskSE(),
52 fCentVar("V0M"),
53 fCentFrom(0),
54 fCentTo(100),
55 fVtxZMin(-10),
56 fVtxZMax(+10),
57 fUseQualFlag(1),
58 fClusName(),
59 fDoNtuple(0),
60 fDoAfterburner(0),
61 fAsymMax(1),
62 fNminCells(1),
63 fMinE(0.100),
64 fMinErat(0),
65 fMinEcc(0),
66 fGeoName("EMCAL_FIRSTYEARV1"),
67 fMinNClusPerTr(50),
68 fIsoDist(0.2),
69 fTrClassNames(""),
70 fTrCuts(0),
71 fPrimTrCuts(0),
e0e1022c 72 fPrimTracksName(""),
4ea96211 73 fDoTrMatGeom(0),
74 fTrainMode(0),
75 fMarkCells(),
76 fMinL0Time(-1),
77 fMaxL0Time(1024),
78 fMcMode(0),
79 fEmbedMode(0),
80 fGeom(0),
81 fReco(0),
469b2bff 82 fTrigName(),
4ea96211 83 fDoPSel(kTRUE),
84 fIsGeoMatsSet(0),
85 fNEvs(0),
86 fOutput(0),
87 fTrClassNamesArr(0),
88 fEsdEv(0),
89 fAodEv(0),
90 fRecPoints(0),
91 fDigits(0),
92 fEsdClusters(0),
93 fEsdCells(0),
94 fAodClusters(0),
95 fAodCells(0),
96 fPtRanges(0),
97 fSelTracks(0),
98 fSelPrimTracks(0),
4ea96211 99 fNtuple(0),
100 fHeader(0),
101 fPrimVert(0),
102 fSpdVert(0),
103 fTpcVert(0),
104 fClusters(0),
105 fTriggers(0),
106 fMcParts(0),
107 fHCuts(0x0),
108 fHVertexZ(0x0),
109 fHVertexZ2(0x0),
110 fHCent(0x0),
111 fHCentQual(0x0),
112 fHTclsBeforeCuts(0x0),
113 fHTclsAfterCuts(0x0),
114 fHColuRow(0x0),
115 fHColuRowE(0x0),
116 fHCellMult(0x0),
117 fHCellE(0x0),
118 fHCellH(0x0),
119 fHCellM(0x0),
120 fHCellM2(0x0),
121 fHCellFreqNoCut(0x0),
122 fHCellFreqCut100M(0x0),
123 fHCellFreqCut300M(0x0),
124 fHCellFreqE(0x0),
125 fHCellCheckE(0x0),
126 fHClustEccentricity(0),
127 fHClustEtaPhi(0x0),
128 fHClustEnergyPt(0x0),
129 fHClustEnergySigma(0x0),
130 fHClustSigmaSigma(0x0),
131 fHClustNCellEnergyRatio(0x0),
132 fHClustEnergyNCell(0x0),
133 fHPrimTrackPt(0x0),
134 fHPrimTrackEta(0x0),
135 fHPrimTrackPhi(0x0),
136 fHMatchDr(0x0),
137 fHMatchDz(0x0),
138 fHMatchEp(0x0),
139 fHPionEtaPhi(0x0),
140 fHPionMggPt(0x0),
141 fHPionMggAsym(0x0),
142 fHPionMggDgg(0x0)
143{
144 // Constructor.
145}
146
ea3fd2d5 147//________________________________________________________________________
148AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
149 : AliAnalysisTaskSE(name),
717fe7de 150 fCentVar("V0M"),
151 fCentFrom(0),
152 fCentTo(100),
76332037 153 fVtxZMin(-10),
154 fVtxZMax(+10),
43cfaa06 155 fUseQualFlag(1),
717fe7de 156 fClusName(),
f5d4ab70 157 fDoNtuple(0),
a49742b5 158 fDoAfterburner(0),
f224d35b 159 fAsymMax(1),
a49742b5 160 fNminCells(1),
3c661da5 161 fMinE(0.100),
f224d35b 162 fMinErat(0),
163 fMinEcc(0),
6bf90832 164 fGeoName("EMCAL_FIRSTYEARV1"),
b6c599fe 165 fMinNClusPerTr(50),
296ea9b4 166 fIsoDist(0.2),
b3ee6797 167 fTrClassNames(""),
0ec74551 168 fTrCuts(0),
b6c599fe 169 fPrimTrCuts(0),
e0e1022c 170 fPrimTracksName(""),
b6c599fe 171 fDoTrMatGeom(0),
3a952328 172 fTrainMode(0),
173 fMarkCells(),
174 fMinL0Time(-1),
175 fMaxL0Time(1024),
38727e64 176 fMcMode(0),
cfd7d5b2 177 fEmbedMode(0),
d595acbb 178 fGeom(0),
296ea9b4 179 fReco(0),
469b2bff 180 fTrigName(),
807016ea 181 fDoPSel(kTRUE),
38727e64 182 fIsGeoMatsSet(0),
183 fNEvs(0),
717fe7de 184 fOutput(0),
b3ee6797 185 fTrClassNamesArr(0),
717fe7de 186 fEsdEv(0),
187 fAodEv(0),
43cfaa06 188 fRecPoints(0),
c2fe5f0e 189 fDigits(0),
717fe7de 190 fEsdClusters(0),
191 fEsdCells(0),
192 fAodClusters(0),
286b47a5 193 fAodCells(0),
fa443410 194 fPtRanges(0),
296ea9b4 195 fSelTracks(0),
b6c599fe 196 fSelPrimTracks(0),
788ca675 197 fNtuple(0),
198 fHeader(0),
199 fPrimVert(0),
200 fSpdVert(0),
201 fTpcVert(0),
202 fClusters(0),
3a952328 203 fTriggers(0),
807016ea 204 fMcParts(0),
fa443410 205 fHCuts(0x0),
206 fHVertexZ(0x0),
76332037 207 fHVertexZ2(0x0),
fa443410 208 fHCent(0x0),
209 fHCentQual(0x0),
b3ee6797 210 fHTclsBeforeCuts(0x0),
211 fHTclsAfterCuts(0x0),
d595acbb 212 fHColuRow(0x0),
213 fHColuRowE(0x0),
214 fHCellMult(0x0),
215 fHCellE(0x0),
216 fHCellH(0x0),
6eb6260e 217 fHCellM(0x0),
a49742b5 218 fHCellM2(0x0),
296ea9b4 219 fHCellFreqNoCut(0x0),
2e4d8148 220 fHCellFreqCut100M(0x0),
221 fHCellFreqCut300M(0x0),
222 fHCellFreqE(0x0),
296ea9b4 223 fHCellCheckE(0x0),
6eb6260e 224 fHClustEccentricity(0),
fa443410 225 fHClustEtaPhi(0x0),
226 fHClustEnergyPt(0x0),
227 fHClustEnergySigma(0x0),
d595acbb 228 fHClustSigmaSigma(0x0),
6eb6260e 229 fHClustNCellEnergyRatio(0x0),
f5e0f1e2 230 fHClustEnergyNCell(0x0),
231 fHPrimTrackPt(0x0),
232 fHPrimTrackEta(0x0),
233 fHPrimTrackPhi(0x0),
b6c599fe 234 fHMatchDr(0x0),
235 fHMatchDz(0x0),
236 fHMatchEp(0x0),
fa443410 237 fHPionEtaPhi(0x0),
238 fHPionMggPt(0x0),
6eb6260e 239 fHPionMggAsym(0x0),
240 fHPionMggDgg(0x0)
ea3fd2d5 241{
717fe7de 242 // Constructor.
ea3fd2d5 243
ea3fd2d5 244 DefineOutput(1, TList::Class());
469b2bff 245 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks,EMCALTrigger.,SPDPileupVertices,TrkPileupVertices "
296ea9b4 246 "AOD:header,vertices,emcalCells,tracks";
ea3fd2d5 247}
ea3fd2d5 248
717fe7de 249//________________________________________________________________________
250AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
251{
252 // Destructor.
ea3fd2d5 253
b3ee6797 254 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
255 delete fOutput; fOutput = 0;
256 }
fa443410 257 delete fPtRanges; fPtRanges = 0;
a2de5ca1 258 fGeom = 0; // do not delete geometry when using instance
296ea9b4 259 delete fReco; fReco = 0;
b3ee6797 260 delete fTrClassNamesArr;
296ea9b4 261 delete fSelTracks;
b6c599fe 262 delete fSelPrimTracks;
d595acbb 263 delete [] fHColuRow;
264 delete [] fHColuRowE;
265 delete [] fHCellMult;
296ea9b4 266 delete [] fHCellFreqNoCut;
2e4d8148 267 delete [] fHCellFreqCut100M;
268 delete [] fHCellFreqCut300M;
ea3fd2d5 269}
717fe7de 270
ea3fd2d5 271//________________________________________________________________________
272void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
273{
717fe7de 274 // Create user objects here.
ea3fd2d5 275
f3582e89 276 cout << "AliAnalysisTaskEMCALPi0PbPb: Input settings" << endl;
277 cout << " fCentVar: " << fCentVar << endl;
278 cout << " fCentFrom: " << fCentFrom << endl;
279 cout << " fCentTo: " << fCentTo << endl;
280 cout << " fVtxZMin: " << fVtxZMin << endl;
281 cout << " fVtxZMax: " << fVtxZMax << endl;
282 cout << " fUseQualFlag: " << fUseQualFlag << endl;
283 cout << " fClusName: \"" << fClusName << "\"" << endl;
284 cout << " fDoNtuple: " << fDoNtuple << endl;
285 cout << " fDoAfterburner: " << fDoAfterburner << endl;
286 cout << " fAsymMax: " << fAsymMax << endl;
287 cout << " fNminCells: " << fNminCells << endl;
288 cout << " fMinE: " << fMinE << endl;
289 cout << " fMinErat: " << fMinErat << endl;
290 cout << " fMinEcc: " << fMinEcc << endl;
291 cout << " fGeoName: \"" << fGeoName << "\"" << endl;
292 cout << " fMinNClusPerTr: " << fMinNClusPerTr << endl;
38727e64 293 cout << " fIsoDist: " << fIsoDist << endl;
f3582e89 294 cout << " fTrClassNames: \"" << fTrClassNames << "\"" << endl;
f3582e89 295 cout << " fTrCuts: " << fTrCuts << endl;
296 cout << " fPrimTrCuts: " << fPrimTrCuts << endl;
297 cout << " fDoTrMatGeom: " << fDoTrMatGeom << endl;
298 cout << " fTrainMode: " << fTrainMode << endl;
299 cout << " fMarkCells: " << fMarkCells << endl;
300 cout << " fMinL0Time: " << fMinL0Time << endl;
301 cout << " fMaxL0Time: " << fMaxL0Time << endl;
38727e64 302 cout << " fMcMode: " << fMcMode << endl;
cfd7d5b2 303 cout << " fEmbedMode: " << fEmbedMode << endl;
38727e64 304 cout << " fGeom: " << fGeom << endl;
305 cout << " fReco: " << fReco << endl;
2ee7bda4 306 cout << " fTrigName: " << fTrigName << endl;
807016ea 307 cout << " fDoPSel: " << fDoPSel << endl;
38727e64 308
309 if (!fGeom)
80d69b70 310 fGeom = AliEMCALGeometry::GetInstance(fGeoName);
38727e64 311 else {
312 if (fGeom->GetMatrixForSuperModule(0))
313 fIsGeoMatsSet = kTRUE;
314 }
315 if (!fReco)
316 fReco = new AliEMCALRecoUtils();
b3ee6797 317 fTrClassNamesArr = fTrClassNames.Tokenize(" ");
717fe7de 318 fOutput = new TList();
4ea96211 319 fOutput->SetOwner(1);
296ea9b4 320 fSelTracks = new TObjArray;
b6c599fe 321 fSelPrimTracks = new TObjArray;
ea3fd2d5 322
f5d4ab70 323 if (fDoNtuple) {
a9c1389d 324 TFile *f = OpenFile(1);
325 TDirectory::TContext context(f);
6bf90832 326 if (f) {
327 f->SetCompressionLevel(2);
788ca675 328 fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
6bf90832 329 fNtuple->SetDirectory(f);
3a952328 330 if (fTrainMode) {
331 fNtuple->SetAutoFlush(-2*1024*1024);
332 fNtuple->SetAutoSave(0);
333 } else {
334 fNtuple->SetAutoFlush(-32*1024*1024);
335 fNtuple->SetAutoSave(0);
336 }
337
788ca675 338 fHeader = new AliStaHeader;
339 fNtuple->Branch("header", &fHeader, 16*1024, 99);
340 fPrimVert = new AliStaVertex;
341 fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
342 fSpdVert = new AliStaVertex;
343 fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
344 fTpcVert = new AliStaVertex;
345 fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
346 if (TClass::GetClass("AliStaCluster"))
347 TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
3a952328 348 fClusters = new TClonesArray("AliStaCluster");
788ca675 349 fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
3a952328 350 if (TClass::GetClass("AliStaTrigger"))
351 TClass::GetClass("AliStaTrigger")->IgnoreTObjectStreamer();
352 fTriggers = new TClonesArray("AliStaTrigger");
353 fNtuple->Branch("l0prim", &fTriggers, 16*1024, 99);
cfd7d5b2 354 if (fMcMode||fEmbedMode) {
807016ea 355 if (TClass::GetClass("AliStaPart"))
356 TClass::GetClass("AliStaPart")->IgnoreTObjectStreamer();
357 fMcParts = new TClonesArray("AliStaPart");
358 fNtuple->Branch("mcparts", &fMcParts, 8*16*1024, 99);
359 }
6bf90832 360 }
f5d4ab70 361 }
362
002dcebe 363 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 364 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
365 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
002dcebe 366
d595acbb 367 // histograms
587ce6c6 368 Bool_t th1 = TH1::GetDefaultSumw2();
a49742b5 369 TH1::SetDefaultSumw2(kTRUE);
587ce6c6 370 Bool_t th2 = TH2::GetDefaultSumw2();
a49742b5 371 TH2::SetDefaultSumw2(kTRUE);
b3ee6797 372 fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
373 fHCuts->GetXaxis()->SetBinLabel(1,"All");
374 fHCuts->GetXaxis()->SetBinLabel(2,"PS");
375 fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
376 fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
377 fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
fa443410 378 fOutput->Add(fHCuts);
379 fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
380 fHVertexZ->SetXTitle("z [cm]");
381 fOutput->Add(fHVertexZ);
76332037 382 fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
383 fHVertexZ2->SetXTitle("z [cm]");
384 fOutput->Add(fHVertexZ2);
f2b8fca6 385 fHCent = new TH1F("hCentBeforeCut","",102,-1,101);
fa443410 386 fHCent->SetXTitle(fCentVar.Data());
76332037 387 fOutput->Add(fHCent);
f2b8fca6 388 fHCentQual = new TH1F("hCentAfterCut","",102,-1,101);
fa443410 389 fHCentQual->SetXTitle(fCentVar.Data());
390 fOutput->Add(fHCentQual);
2e4d8148 391 fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
392 fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
b3ee6797 393 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
394 const char *name = fTrClassNamesArr->At(i)->GetName();
395 fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
396 fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
397 }
398 fOutput->Add(fHTclsBeforeCuts);
399 fOutput->Add(fHTclsAfterCuts);
90d5b88b 400
d595acbb 401 // histograms for cells
296ea9b4 402 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
403 fHColuRow = new TH2*[nsm];
404 fHColuRowE = new TH2*[nsm];
405 fHCellMult = new TH1*[nsm];
406 for (Int_t i = 0; i<nsm; ++i) {
2e4d8148 407 fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
d595acbb 408 fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
409 fHColuRow[i]->SetXTitle("col (i#eta)");
410 fHColuRow[i]->SetYTitle("row (i#phi)");
2e4d8148 411 fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
90d5b88b 412 fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
d595acbb 413 fHColuRowE[i]->SetXTitle("col (i#eta)");
414 fHColuRowE[i]->SetYTitle("row (i#phi)");
415 fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
416 fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
90d5b88b 417 fHCellMult[i]->SetXTitle("# of cells");
d595acbb 418 fOutput->Add(fHColuRow[i]);
419 fOutput->Add(fHColuRowE[i]);
420 fOutput->Add(fHCellMult[i]);
421 }
2e4d8148 422 fHCellE = new TH1F("hCellE","",250,0.,25.);
d595acbb 423 fHCellE->SetXTitle("E_{cell} [GeV]");
424 fOutput->Add(fHCellE);
fb5a0b69 425 fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
6eb6260e 426 fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
d595acbb 427 fOutput->Add(fHCellH);
fb5a0b69 428 fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
6eb6260e 429 fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
430 fOutput->Add(fHCellM);
fb5a0b69 431 fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
f4ec884e 432 fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
a49742b5 433 fOutput->Add(fHCellM2);
90d5b88b 434
2e4d8148 435 fHCellFreqNoCut = new TH1*[nsm];
436 fHCellFreqCut100M = new TH1*[nsm];
437 fHCellFreqCut300M = new TH1*[nsm];
438 fHCellFreqE = new TH1*[nsm];
296ea9b4 439 for (Int_t i = 0; i<nsm; ++i){
440 Double_t lbin = i*24*48-0.5;
441 Double_t hbin = lbin+24*48;
2e4d8148 442 fHCellFreqNoCut[i] = new TH1F(Form("hCellFreqNoCut_SM%d",i),
443 Form("Frequency SM%d (no cut);id;#",i), 1152, lbin, hbin);
444 fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
445 Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
446 fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
447 Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
448 fHCellFreqE[i] = new TH1F(Form("hCellFreqE_SM%d",i),
449 Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
296ea9b4 450 fOutput->Add(fHCellFreqNoCut[i]);
2e4d8148 451 fOutput->Add(fHCellFreqCut100M[i]);
452 fOutput->Add(fHCellFreqCut300M[i]);
453 fOutput->Add(fHCellFreqE[i]);
296ea9b4 454 }
3a952328 455 if (!fMarkCells.IsNull()) {
296ea9b4 456 fHCellCheckE = new TH1*[24*48*nsm];
408dc04c 457 memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
3a952328 458 TObjArray *cells = fMarkCells.Tokenize(" ");
459 Int_t n = cells->GetEntries();
460 Int_t *tcs = new Int_t[n];
461 for (Int_t i=0;i<n;++i) {
462 TString name(cells->At(i)->GetName());
463 tcs[i]=name.Atoi();
464 }
465 for (Int_t i = 0; i<n; ++i) {
296ea9b4 466 Int_t c=tcs[i];
467 if (c<24*48*nsm) {
3a952328 468 fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 1000, 0, 10);
296ea9b4 469 fOutput->Add(fHCellCheckE[i]);
470 }
471 }
3a952328 472 delete cells;
473 delete [] tcs;
296ea9b4 474 }
475
d595acbb 476 // histograms for clusters
3a952328 477 if (!fTrainMode) {
478 fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
479 fHClustEccentricity->SetXTitle("#epsilon_{C}");
480 fOutput->Add(fHClustEccentricity);
481 fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
482 fHClustEtaPhi->SetXTitle("#eta");
483 fHClustEtaPhi->SetYTitle("#varphi");
484 fOutput->Add(fHClustEtaPhi);
485 fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
486 fHClustEnergyPt->SetXTitle("E [GeV]");
487 fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
488 fOutput->Add(fHClustEnergyPt);
489 fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
490 fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
491 fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
492 fOutput->Add(fHClustEnergySigma);
493 fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
494 fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
495 fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
496 fOutput->Add(fHClustSigmaSigma);
497 fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
498 fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
499 fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
500 fOutput->Add(fHClustNCellEnergyRatio);
f5e0f1e2 501 fHClustEnergyNCell = new TH2F("hClustEnergyNCell","",200,0,100,50,0,50);
502 fHClustEnergyNCell->SetXTitle("E_{clus}");
503 fHClustEnergyNCell->SetYTitle("N_{cells}");
504 fOutput->Add(fHClustEnergyNCell);
3a952328 505 }
90d5b88b 506
f5e0f1e2 507 // histograms for primary tracks
508 fHPrimTrackPt = new TH1F("hPrimTrackPt",";p_{T} [GeV/c]",500,0,50);
509 fOutput->Add(fHPrimTrackPt);
510 fHPrimTrackEta = new TH1F("hPrimTrackEta",";#eta",40,-2,2);
511 fOutput->Add(fHPrimTrackEta);
512 fHPrimTrackPhi = new TH1F("hPrimTrackPhi",";#varPhi [rad]",63,0,6.3);
513 fOutput->Add(fHPrimTrackPhi);
b6c599fe 514 // histograms for track matching
515 fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
516 fOutput->Add(fHMatchDr);
517 fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
518 fOutput->Add(fHMatchDz);
519 fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
520 fOutput->Add(fHMatchEp);
3a952328 521
d595acbb 522 // histograms for pion candidates
3a952328 523 if (!fTrainMode) {
524 fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
525 fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
526 fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
527 fOutput->Add(fHPionEtaPhi);
528 fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
529 fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
530 fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
531 fOutput->Add(fHPionMggPt);
532 fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
533 fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
534 fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
535 fOutput->Add(fHPionMggAsym);
536 fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
537 fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
538 fHPionMggDgg->SetYTitle("opening angle [grad]");
539 fOutput->Add(fHPionMggDgg);
540 const Int_t nbins = 20;
541 Double_t xbins[nbins] = {0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,9,10,12.5,15,20,25,50};
542 fPtRanges = new TAxis(nbins-1,xbins);
543 for (Int_t i = 0; i<=nbins; ++i) {
544 fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
545 fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
546 if (i==0)
547 fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
548 else if (i==nbins)
549 fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
550 else
551 fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
552 fOutput->Add(fHPionInvMasses[i]);
553 }
fa443410 554 }
296ea9b4 555
587ce6c6 556 TH1::SetDefaultSumw2(th1);
557 TH2::SetDefaultSumw2(th2);
717fe7de 558 PostData(1, fOutput);
ea3fd2d5 559}
560
561//________________________________________________________________________
562void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
563{
717fe7de 564 // Called for each event.
565
566 if (!InputEvent())
567 return;
568
43cfaa06 569 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
570 fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
44310bce 571 UInt_t offtrigger = 0;
43cfaa06 572 if (fEsdEv) {
573 am->LoadBranch("AliESDRun.");
574 am->LoadBranch("AliESDHeader.");
9809ed8c 575 UInt_t mask1 = fEsdEv->GetESDRun()->GetDetectorsInDAQ();
576 UInt_t mask2 = fEsdEv->GetESDRun()->GetDetectorsInReco();
577 Bool_t desc1 = (mask1 >> 18) & 0x1;
578 Bool_t desc2 = (mask2 >> 18) & 0x1;
2cee9afd 579 if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(18)=="EMCAL"
9809ed8c 580 AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
581 mask1, fEsdEv->GetESDRun()->GetDetectorsInReco(),
582 mask2, fEsdEv->GetESDRun()->GetDetectorsInDAQ()));
44310bce 583 return;
584 }
585 offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
43cfaa06 586 } else {
587 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
408dc04c 588 if (!fAodEv) {
589 AliFatal("Neither ESD nor AOD event found");
590 return;
591 }
43cfaa06 592 am->LoadBranch("header");
44310bce 593 offtrigger = fAodEv->GetHeader()->GetOfflineTrigger();
594 }
38727e64 595 if (!fMcMode && (offtrigger & AliVEvent::kFastOnly)) {
44310bce 596 AliWarning(Form("EMCAL not in fast only partition"));
597 return;
43cfaa06 598 }
b6c599fe 599 if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
27c2e3d9 600 AliWarning("Accessing geometry from OCDB, this is not very efficient!");
601 AliCDBManager *cdb = AliCDBManager::Instance();
602 if (!cdb->IsDefaultStorageSet())
603 cdb->SetDefaultStorage("raw://");
604 Int_t runno = InputEvent()->GetRunNumber();
605 if (runno != cdb->GetRun())
606 cdb->SetRun(runno);
607 AliGeomManager::LoadGeometry();
608 }
609
610 if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
611 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
9809ed8c 612 for (Int_t i=0; i<nsm; ++i) {
613 const TGeoHMatrix *geom = 0;
614 if (fEsdEv)
615 geom = fEsdEv->GetESDRun()->GetEMCALMatrix(i);
616 else
617 geom = fAodEv->GetHeader()->GetEMCALMatrix(i);
618 if (!geom)
619 continue;
620 geom->Print();
621 fGeom->SetMisalMatrix(geom,i);
27c2e3d9 622 }
623 fIsGeoMatsSet = kTRUE;
624 }
625
626 if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
627 if (fEsdEv) {
628 const AliESDRun *erun = fEsdEv->GetESDRun();
629 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
630 erun->GetCurrentDip(),
631 AliMagF::kConvLHC,
632 kFALSE,
633 erun->GetBeamEnergy(),
634 erun->GetBeamType());
635 TGeoGlobalMagField::Instance()->SetField(field);
636 } else {
637 Double_t pol = -1; //polarity
638 Double_t be = -1; //beam energy
639 AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
640 Int_t runno = fAodEv->GetRunNumber();
641 if (runno>=136851 && runno<138275) {
642 pol = -1;
643 be = 2760;
644 btype = AliMagF::kBeamTypeAA;
645 } else if (runno>=138275 && runno<=139517) {
646 pol = +1;
647 be = 2760;
648 btype = AliMagF::kBeamTypeAA;
649 } else {
650 AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
651 }
652 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
653 }
654 }
655
b3ee6797 656 Int_t cut = 1;
657 fHCuts->Fill(cut++);
658
659 TString trgclasses;
660 AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
661 if (h) {
662 trgclasses = fEsdEv->GetFiredTriggerClasses();
663 } else {
664 AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
665 if (h2)
666 trgclasses = h2->GetFiredTriggerClasses();
667 }
668 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
669 const char *name = fTrClassNamesArr->At(i)->GetName();
0ceca1f8 670 TRegexp regexp(name);
671 if (trgclasses.Contains(regexp))
b3ee6797 672 fHTclsBeforeCuts->Fill(1+i);
673 }
674
807016ea 675 if (fDoPSel && offtrigger==0)
b3ee6797 676 return;
677
fa443410 678 fHCuts->Fill(cut++);
286b47a5 679
717fe7de 680 const AliCentrality *centP = InputEvent()->GetCentrality();
681 Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
fa443410 682 fHCent->Fill(cent);
717fe7de 683 if (cent<fCentFrom||cent>fCentTo)
286b47a5 684 return;
43cfaa06 685
fa443410 686 fHCuts->Fill(cut++);
286b47a5 687
43cfaa06 688 if (fUseQualFlag) {
689 if (centP->GetQuality()>0)
690 return;
691 }
286b47a5 692
fa443410 693 fHCentQual->Fill(cent);
694 fHCuts->Fill(cut++);
717fe7de 695
43cfaa06 696 if (fEsdEv) {
fa443410 697 am->LoadBranch("PrimaryVertex.");
698 am->LoadBranch("SPDVertex.");
699 am->LoadBranch("TPCVertex.");
43cfaa06 700 } else {
701 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
702 am->LoadBranch("vertices");
3c661da5 703 if (!fAodEv) return;
43cfaa06 704 }
705
706 const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
717fe7de 707 if (!vertex)
708 return;
709
fa443410 710 fHVertexZ->Fill(vertex->GetZ());
286b47a5 711
717fe7de 712 if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
713 return;
286b47a5 714
fa443410 715 fHCuts->Fill(cut++);
76332037 716 fHVertexZ2->Fill(vertex->GetZ());
717fe7de 717
b3ee6797 718 // count number of accepted events
719 ++fNEvs;
720
721 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
722 const char *name = fTrClassNamesArr->At(i)->GetName();
0ceca1f8 723 TRegexp regexp(name);
724 if (trgclasses.Contains(regexp))
b3ee6797 725 fHTclsAfterCuts->Fill(1+i);
726 }
727
43cfaa06 728 fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
c2fe5f0e 729 fDigits = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
b6c599fe 730 fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 731 fEsdCells = 0; // will be set if ESD input used
b6c599fe 732 fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
43cfaa06 733 // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
734 fAodCells = 0; // will be set if AOD input used
717fe7de 735
43cfaa06 736 // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
2ee7bda4 737 Bool_t overwrite = 0;
43cfaa06 738 Bool_t clusattached = 0;
739 Bool_t recalibrated = 0;
740 if (1 && !fClusName.IsNull()) {
741 AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
742 TObjArray *ts = am->GetTasks();
743 cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
744 if (cltask && cltask->GetClusters()) {
2ee7bda4 745 fRecPoints = cltask->GetClusters();
746 fDigits = cltask->GetDigits();
43cfaa06 747 clusattached = cltask->GetAttachClusters();
2ee7bda4 748 overwrite = cltask->GetOverwrite();
43cfaa06 749 if (cltask->GetCalibData()!=0)
750 recalibrated = kTRUE;
751 }
752 }
38727e64 753 if (1 && !fClusName.IsNull()) {
754 TList *l = 0;
755 if (AODEvent())
756 l = AODEvent()->GetList();
757 else if (fAodEv)
758 l = fAodEv->GetList();
717fe7de 759 if (l) {
38727e64 760 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
717fe7de 761 }
ea3fd2d5 762 }
717fe7de 763
764 if (fEsdEv) { // ESD input mode
43cfaa06 765 if (1 && (!fRecPoints||clusattached)) {
2ee7bda4 766 if (!clusattached && !overwrite)
43cfaa06 767 am->LoadBranch("CaloClusters");
717fe7de 768 TList *l = fEsdEv->GetList();
2ee7bda4 769 if (clusattached) {
770 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
1829ebfe 771 } else {
38727e64 772 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
ea3fd2d5 773 }
774 }
43cfaa06 775 if (1) {
776 if (!recalibrated)
777 am->LoadBranch("EMCALCells.");
717fe7de 778 fEsdCells = fEsdEv->GetEMCALCells();
779 }
780 } else if (fAodEv) { // AOD input mode
43cfaa06 781 if (1 && (!fAodClusters || clusattached)) {
782 if (!clusattached)
783 am->LoadBranch("caloClusters");
717fe7de 784 TList *l = fAodEv->GetList();
717fe7de 785 if (l) {
38727e64 786 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
ea3fd2d5 787 }
788 }
43cfaa06 789 if (1) {
790 if (!recalibrated)
791 am->LoadBranch("emcalCells");
717fe7de 792 fAodCells = fAodEv->GetEMCALCells();
793 }
794 } else {
795 AliFatal("Impossible to not have either pointer to ESD or AOD event");
796 }
ea3fd2d5 797
43cfaa06 798 if (1) {
799 AliDebug(2,Form("fRecPoints set: %p", fRecPoints));
c2fe5f0e 800 AliDebug(2,Form("fDigits set: %p", fDigits));
43cfaa06 801 AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
802 AliDebug(2,Form("fEsdCells set: %p", fEsdCells));
803 AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
804 AliDebug(2,Form("fAodCells set: %p", fAodCells));
805 }
806
a49742b5 807 if (fDoAfterburner)
808 ClusterAfterburner();
6eb6260e 809
38727e64 810 CalcMcInfo();
3a952328 811 CalcCaloTriggers();
b6c599fe 812 CalcPrimTracks();
3a952328 813 CalcTracks();
296ea9b4 814 CalcClusterProps();
815
286b47a5 816 FillCellHists();
3a952328 817 if (!fTrainMode) {
818 FillClusHists();
819 FillPionHists();
820 FillOtherHists();
821 }
38727e64 822 FillMcHists();
788ca675 823 FillNtuple();
ea3fd2d5 824
3a952328 825 if (fTrainMode) {
826 fSelTracks->Clear();
827 fSelPrimTracks->Clear();
cfd7d5b2 828 if (fMcParts)
807016ea 829 fMcParts->Clear();
4ea96211 830 if (fTriggers)
831 fTriggers->Clear();
832 if (fClusters)
833 fClusters->Clear();
3a952328 834 }
de4cee41 835
717fe7de 836 PostData(1, fOutput);
ea3fd2d5 837}
838
839//________________________________________________________________________
840void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
841{
717fe7de 842 // Terminate called at the end of analysis.
f5d4ab70 843
844 if (fNtuple) {
4ea96211 845 TFile *f = OpenFile(1);
a9c1389d 846 TDirectory::TContext context(f);
f5d4ab70 847 if (f)
848 fNtuple->Write();
849 }
d9f26424 850
bc8a45bd 851 AliInfo(Form("%s: Accepted %lld events ", GetName(), fNEvs));
ea3fd2d5 852}
ea3fd2d5 853
296ea9b4 854//________________________________________________________________________
3a952328 855void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
296ea9b4 856{
3a952328 857 // Calculate triggers
296ea9b4 858
38727e64 859 if (fAodEv)
860 return; // information not available in AOD
861
4ea96211 862 if (!fTriggers)
863 return;
864
3a952328 865 fTriggers->Clear();
296ea9b4 866
2ee7bda4 867 if (fTrigName.Length()<=0)
3a952328 868 return;
869
2ee7bda4 870 TClonesArray *arr = dynamic_cast<TClonesArray*>(fEsdEv->FindListObject(fTrigName));
871 if (!arr) {
872 AliError(Form("Could not get array with name %s", fTrigName.Data()));
3a952328 873 return;
3a952328 874 }
875
2ee7bda4 876 Int_t nNumberOfCaloClusters = arr->GetEntries();
877 for(Int_t j = 0, ntrigs = 0; j < nNumberOfCaloClusters; ++j) {
878 AliVCluster *cl = dynamic_cast<AliVCluster*>(arr->At(j));
879 if (!cl)
3a952328 880 continue;
2ee7bda4 881 if (!cl->IsEMCAL())
3a952328 882 continue;
2ee7bda4 883 if (cl->E()<1)
c47674e0 884 continue;
3a952328 885 AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
2ee7bda4 886 Float_t pos[3] = {0,0,0};
887 cl->GetPosition(pos);
888 TVector3 vpos(pos);
889 trignew->fE = cl->E();
890 trignew->fEta = vpos.Eta();
891 trignew->fPhi = vpos.Phi();
892 Short_t id = -1;
893 GetMaxCellEnergy(cl, id);
894 trignew->fIdMax = id;
3a952328 895 }
896}
897
898//________________________________________________________________________
899void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
900{
901 // Calculate cluster properties
902
4ea96211 903 if (!fClusters)
904 return;
905
3a952328 906 fClusters->Clear();
907
908 AliVCaloCells *cells = fEsdCells;
909 if (!cells)
910 cells = fAodCells;
911 if (!cells)
912 return;
913
914 TObjArray *clusters = fEsdClusters;
915 if (!clusters)
916 clusters = fAodClusters;
917 if (!clusters)
918 return;
919
920 Int_t ncells = cells->GetNumberOfCells();
921 Int_t nclus = clusters->GetEntries();
922 Int_t ntrks = fSelTracks->GetEntries();
3a952328 923
924 std::map<Short_t,Short_t> map;
925 for (Short_t pos=0;pos<ncells;++pos) {
926 Short_t id = cells->GetCellNumber(pos);
927 map[id]=pos;
928 }
929
8c56d760 930 TObjArray filtMcParts;
cfd7d5b2 931 if (fMcParts) {
8c56d760 932 Int_t nmc = fMcParts->GetEntries();
933 for (Int_t i=0; i<nmc; ++i) {
934 AliStaPart *pa = static_cast<AliStaPart*>(fMcParts->At(i));
935 if (pa->OnEmcal())
936 filtMcParts.Add(pa);
937 }
938 }
939
fe38566a 940 Double_t vertex[3] = {0,0,0};
941 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
942
3a952328 943 for(Int_t i=0, ncl=0; i<nclus; ++i) {
944 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
b5ba4932 945
3a952328 946 if (!clus)
947 continue;
948 if (!clus->IsEMCAL())
949 continue;
950 if (clus->E()<fMinE)
951 continue;
952
fe38566a 953 Float_t clsPos[3] = {0,0,0};
3a952328 954 clus->GetPosition(clsPos);
955 TVector3 clsVec(clsPos);
3a952328 956 TLorentzVector clusterVec;
957 clus->GetMomentum(clusterVec,vertex);
958 Double_t clsEta = clusterVec.Eta();
959
960 AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
961 cl->fE = clus->E();
962 cl->fR = clsVec.Perp();
963 cl->fEta = clsVec.Eta();
964 cl->fPhi = clsVec.Phi();
965 cl->fN = clus->GetNCells();
966 cl->fN1 = GetNCells(clus,0.100);
967 cl->fN3 = GetNCells(clus,0.300);
968 Short_t id = -1;
969 Double_t emax = GetMaxCellEnergy(clus, id);
970 cl->fIdMax = id;
2ee7bda4 971 cl->fSM = fGeom->GetSuperModuleNumber(id);
3a952328 972 cl->fEmax = emax;
1f41ed3d 973 Short_t id2 = -1;
974 cl->fE2max = GetSecondMaxCellEnergy(clus,id2);
2ee7bda4 975 cl->fTmax = cells->GetCellTime(id);
3a952328 976 if (clus->GetDistanceToBadChannel()<10000)
977 cl->fDbc = clus->GetDistanceToBadChannel();
978 if (!TMath::IsNaN(clus->GetDispersion()))
979 cl->fDisp = clus->GetDispersion();
980 if (!TMath::IsNaN(clus->GetM20()))
0fbe8d4f 981 cl->fM20 = clus->GetM20();
3a952328 982 if (!TMath::IsNaN(clus->GetM02()))
0fbe8d4f 983 cl->fM02 = clus->GetM02();
2ee7bda4 984 Double_t maxAxis = -1, minAxis = -1;
3a952328 985 GetSigma(clus,maxAxis,minAxis);
2ee7bda4 986 clus->SetTOF(maxAxis); // store sigma in TOF for later plotting
3a952328 987 cl->fSig = maxAxis;
2ee7bda4 988 Double_t sEtaEta = -1;
989 Double_t sPhiPhi = -1;
5fc7508c 990 GetSigmaEtaEta(clus, sEtaEta, sPhiPhi);
991 cl->fSigEtaEta = sEtaEta;
992 cl->fSigPhiPhi = sPhiPhi;
2ee7bda4 993 Double_t clusterEcc = -1;
3a952328 994 if (maxAxis > 0)
995 clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
2ee7bda4 996 clus->SetChi2(clusterEcc); // store ecc in chi2 for later plotting
997 cl->fEcc = clusterEcc;
998 cl->fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
999 cl->fTrIso1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
1000 cl->fTrIso2 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
5fc7508c 1001 cl->fTrIsoD1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1);
1002 cl->fTrIso1D1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 1);
1003 cl->fTrIso2D1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist-0.1, 2);
1004 cl->fTrIsoD3 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1);
1005 cl->fTrIso1D3 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 1);
1006 cl->fTrIso2D3 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.1, 2);
1f41ed3d 1007 cl->fTrIsoD4 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2);
1008 cl->fTrIso1D4 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2, 1);
1009 cl->fTrIso2D4 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist+0.2, 2);
5fc7508c 1010 cl->fTrIsoStrip = GetTrackIsoStrip(clusterVec.Eta(), clusterVec.Phi());
2ee7bda4 1011 cl->fCeCore = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
1012 cl->fCeIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
1013 cl->fCeIso1 = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.10);
1014 cl->fCeIso3 = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.30);
1f41ed3d 1015 cl->fCeIso4 = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.40);
2ee7bda4 1016 cl->fCeIso4x4 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 4, 4);
1017 cl->fCeIso5x5 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 5, 5);
1018 cl->fCeIso3x22 = GetCellIsoNxM(clsVec.Eta(),clsVec.Phi(), 3, 22);
1019 cl->fIsShared = IsShared(clus);
1020 cl->fTrigId = -1;
1021 cl->fTrigE = 0;
e0e1022c 1022
2ee7bda4 1023 if (fTriggers) {
1024 Int_t ntrig = fTriggers->GetEntries();
1025 for (Int_t j = 0; j<ntrig; ++j) {
1026 AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1027 if (!sta)
38727e64 1028 continue;
2ee7bda4 1029 Short_t idmax = sta->fIdMax;
1030 Bool_t inc = IsIdPartOfCluster(clus, idmax);
1031 if (inc) {
1032 cl->fTrigId = j;
1033 cl->fTrigE = sta->fE;
1034 break;
1035 }
38727e64 1036 }
3a952328 1037 }
1038
1039 // track matching
1040 Double_t mind2 = 1e10;
1041 for(Int_t j = 0; j<ntrks; ++j) {
1042 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
b6c599fe 1043 if (!track)
1044 continue;
3a952328 1045
1046 if (TMath::Abs(clsEta-track->Eta())>0.5)
b6c599fe 1047 continue;
3a952328 1048
1049 TVector3 vec(clsPos);
b6c599fe 1050
e0e1022c 1051 Float_t tmpEta=-1, tmpPhi=-1, tmpR2=1e10;
1052 Double_t dedx = -1;
1053 if (fEsdEv) {
1054 AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
1055 dedx = esdTrack->GetTPCsignal();
1056 }
1057 if (fDoTrMatGeom) {
3a952328 1058 if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
1059 continue;
e0e1022c 1060 AliExternalTrackParam tParam;
1061 tParam.CopyFromVTrack(track);
1062 if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpEta, tmpPhi))
3a952328 1063 continue;
e0e1022c 1064 tmpR2 = tmpEta*tmpEta + tmpPhi*tmpPhi;
b6c599fe 1065 }
e0e1022c 1066 Double_t d2 = tmpR2;
3a952328 1067 if (mind2>d2) {
1068 mind2=d2;
e0e1022c 1069 cl->fTrDz = tmpEta;
1070 cl->fTrDr = tmpPhi;
3a952328 1071 cl->fTrEp = clus->E()/track->P();
5fc7508c 1072 cl->fTrDedx = dedx;
f3582e89 1073 cl->fIsTrackM = 1;
3a952328 1074 }
1075 }
1076
f3582e89 1077 if (cl->fIsTrackM) {
3a952328 1078 fHMatchDr->Fill(cl->fTrDr);
1079 fHMatchDz->Fill(cl->fTrDz);
1080 fHMatchEp->Fill(cl->fTrEp);
b6c599fe 1081 }
8c56d760 1082
1083 //mc matching
cfd7d5b2 1084 if (fMcParts) {
8c56d760 1085 Int_t nmc = filtMcParts.GetEntries();
1086 Double_t diffR2 = 1e9;
1087 AliStaPart *msta = 0;
1088 for (Int_t j=0; j<nmc; ++j) {
1089 AliStaPart *pa = static_cast<AliStaPart*>(filtMcParts.At(j));
1090 Double_t t1=clsVec.Eta()-pa->fVEta;
1091 Double_t t2=TVector2::Phi_mpi_pi(clsVec.Phi()-pa->fVPhi);
1092 Double_t tmp = t1*t1+t2*t2;
1093 if (tmp<diffR2) {
1094 diffR2 = tmp;
1095 msta = pa;
1096 }
1097 }
c2fe5f0e 1098 if (diffR2<10 && msta!=0) {
1099 cl->fMcLabel = msta->fLab;
1100 }
1101 }
1102
1103 cl->fEmbE = 0;
1104 if (fDigits && fEmbedMode) {
1105 for(Int_t j=0; j<cl->fN; ++j) {
1106 Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
1107 Short_t pos = -1;
1108 std::map<Short_t,Short_t>::iterator it = map.find(cid);
1109 if (it!=map.end())
1110 pos = it->second;
1111 if (pos<0)
1112 continue;
1113 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
1114 if (!digit)
1115 continue;
1116 if (digit->GetId() != cid) {
1117 AliError(Form("Ids should be equal: %d %d", cid, digit->GetId()));
1118 continue;
1119 }
1120 if (digit->GetType()<-1) {
1121 cl->fEmbE += digit->GetChi2();
1122 }
1123 }
8c56d760 1124 }
b6c599fe 1125 }
1126}
1127
1128//________________________________________________________________________
1129void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
1130{
a2de5ca1 1131 // Calculate track properties for primary tracks.
b6c599fe 1132
1133 fSelPrimTracks->Clear();
1134
1135 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1136 TClonesArray *tracks = 0;
e0e1022c 1137 Bool_t doConstrain = kFALSE;
1138 TString trkname(fPrimTracksName);
b6c599fe 1139 if (fEsdEv) {
e0e1022c 1140 if (trkname.IsNull()) {
1141 trkname = "Tracks";
1142 am->LoadBranch("Tracks");
1143 fSelPrimTracks->SetOwner(kTRUE);
1144 doConstrain = kTRUE;
1145 }
b6c599fe 1146 TList *l = fEsdEv->GetList();
e0e1022c 1147 tracks = dynamic_cast<TClonesArray*>(l->FindObject(trkname));
b6c599fe 1148 } else {
e0e1022c 1149 trkname = "tracks";
b6c599fe 1150 am->LoadBranch("tracks");
1151 TList *l = fAodEv->GetList();
1152 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1153 }
1154
e0e1022c 1155 if (!tracks) {
1156 AliError(Form("Pointer to tracks %s == 0", trkname.Data() ));
296ea9b4 1157 return;
e0e1022c 1158 }
296ea9b4 1159
0ec74551 1160 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
b6c599fe 1161 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
1162 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
0ec74551 1163
296ea9b4 1164 if (fEsdEv) {
0ec74551 1165 am->LoadBranch("PrimaryVertex.");
1166 const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
1167 am->LoadBranch("SPDVertex.");
1168 const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
1169 am->LoadBranch("Tracks");
1170 const Int_t Ntracks = tracks->GetEntries();
1171 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1172 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1173 if (!track) {
1174 AliWarning(Form("Could not receive track %d\n", iTracks));
1175 continue;
1176 }
e0e1022c 1177 if (fPrimTrCuts && !fPrimTrCuts->IsSelected(track))
296ea9b4 1178 continue;
1179 Double_t eta = track->Eta();
1180 if (eta<-1||eta>1)
1181 continue;
0ec74551 1182 if (track->Phi()<phimin||track->Phi()>phimax)
1183 continue;
2e4d8148 1184
e0e1022c 1185 if (!doConstrain) {
1186 fSelPrimTracks->Add(track);
1187 continue;
1188 }
0ec74551 1189 AliESDtrack copyt(*track);
1190 Double_t bfield[3];
1191 copyt.GetBxByBz(bfield);
1192 AliExternalTrackParam tParam;
1193 Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1194 if (!relate)
1195 continue;
1196 copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
1197 Double_t p[3] = { 0. };
1198 copyt.GetPxPyPz(p);
1199 Double_t pos[3] = { 0. };
1200 copyt.GetXYZ(pos);
1201 Double_t covTr[21] = { 0. };
1202 copyt.GetCovarianceXYZPxPyPz(covTr);
1203 Double_t pid[10] = { 0. };
1204 copyt.GetESDpid(pid);
1205 AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1206 copyt.GetLabel(),
1207 p,
1208 kTRUE,
1209 pos,
1210 kFALSE,
1211 covTr,
1212 (Short_t)copyt.GetSign(),
1213 copyt.GetITSClusterMap(),
1214 pid,
1215 0,/*fPrimaryVertex,*/
1216 kTRUE, // check if this is right
1217 vtx->UsesTrack(copyt.GetID()));
1218 aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1219 aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1220 Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1221 if(ndf>0)
1222 aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1223 else
1224 aTrack->SetChi2perNDF(-1);
1225 aTrack->SetFlags(copyt.GetStatus());
1226 aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
b6c599fe 1227 fSelPrimTracks->Add(aTrack);
296ea9b4 1228 }
1229 } else {
1230 Int_t ntracks = tracks->GetEntries();
1231 for (Int_t i=0; i<ntracks; ++i) {
1232 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1233 if (!track)
1234 continue;
1235 Double_t eta = track->Eta();
1236 if (eta<-1||eta>1)
1237 continue;
0ec74551 1238 if (track->Phi()<phimin||track->Phi()>phimax)
1239 continue;
b6c599fe 1240 if(track->GetTPCNcls()<fMinNClusPerTr)
296ea9b4 1241 continue;
b6c599fe 1242 //todo: Learn how to set/filter AODs for prim/sec tracks
1243 fSelPrimTracks->Add(track);
296ea9b4 1244 }
1245 }
1246}
1247
1248//________________________________________________________________________
3a952328 1249void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
296ea9b4 1250{
3a952328 1251 // Calculate track properties (including secondaries).
b6c599fe 1252
3a952328 1253 fSelTracks->Clear();
296ea9b4 1254
3a952328 1255 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1256 TClonesArray *tracks = 0;
1257 if (fEsdEv) {
1258 am->LoadBranch("Tracks");
1259 TList *l = fEsdEv->GetList();
1260 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1261 } else {
1262 am->LoadBranch("tracks");
1263 TList *l = fAodEv->GetList();
1264 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1265 }
296ea9b4 1266
3a952328 1267 if (!tracks)
1268 return;
296ea9b4 1269
3a952328 1270 if (fEsdEv) {
1271 const Int_t Ntracks = tracks->GetEntries();
1272 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1273 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1274 if (!track) {
1275 AliWarning(Form("Could not receive track %d\n", iTracks));
1276 continue;
1277 }
1278 if (fTrCuts && !fTrCuts->IsSelected(track))
1279 continue;
1280 Double_t eta = track->Eta();
1281 if (eta<-1||eta>1)
1282 continue;
1283 fSelTracks->Add(track);
1284 }
1285 } else {
1286 Int_t ntracks = tracks->GetEntries();
1287 for (Int_t i=0; i<ntracks; ++i) {
1288 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
296ea9b4 1289 if (!track)
1290 continue;
3a952328 1291 Double_t eta = track->Eta();
1292 if (eta<-1||eta>1)
c4236a58 1293 continue;
3a952328 1294 if(track->GetTPCNcls()<fMinNClusPerTr)
b6c599fe 1295 continue;
2e4d8148 1296
3a952328 1297 if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
e0e1022c 1298 AliExternalTrackParam tParam;
1299 tParam.CopyFromVTrack(track);
3a952328 1300 if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
1301 Double_t trkPos[3];
1302 tParam.GetXYZ(trkPos);
1303 track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1304 }
b6c599fe 1305 }
3a952328 1306 fSelTracks->Add(track);
c4236a58 1307 }
296ea9b4 1308 }
1309}
1310
323834f0 1311//________________________________________________________________________
3a952328 1312void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
323834f0 1313{
296ea9b4 1314 // Run custer reconstruction afterburner.
323834f0 1315
1316 AliVCaloCells *cells = fEsdCells;
1317 if (!cells)
1318 cells = fAodCells;
1319
1320 if (!cells)
1321 return;
1322
1323 Int_t ncells = cells->GetNumberOfCells();
1324 if (ncells<=0)
1325 return;
1326
1327 Double_t cellMeanE = 0, cellSigE = 0;
1328 for (Int_t i = 0; i<ncells; ++i) {
1329 Double_t cellE = cells->GetAmplitude(i);
1330 cellMeanE += cellE;
1331 cellSigE += cellE*cellE;
1332 }
1333 cellMeanE /= ncells;
1334 cellSigE /= ncells;
1335 cellSigE -= cellMeanE*cellMeanE;
1336 if (cellSigE<0)
1337 cellSigE = 0;
1338 cellSigE = TMath::Sqrt(cellSigE / ncells);
1339
1340 Double_t subE = cellMeanE - 7*cellSigE;
1341 if (subE<0)
1342 return;
1343
1344 for (Int_t i = 0; i<ncells; ++i) {
1345 Short_t id=-1;
1346 Double_t amp=0,time=0;
1347 if (!cells->GetCell(i, id, amp, time))
1348 continue;
1349 amp -= cellMeanE;
1350 if (amp<0.001)
1351 amp = 0;
1352 cells->SetCell(i, id, amp, time);
1353 }
1354
1355 TObjArray *clusters = fEsdClusters;
1356 if (!clusters)
1357 clusters = fAodClusters;
323834f0 1358 if (!clusters)
1359 return;
1360
1361 Int_t nclus = clusters->GetEntries();
1362 for (Int_t i = 0; i<nclus; ++i) {
1363 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1364 if (!clus->IsEMCAL())
1365 continue;
1366 Int_t nc = clus->GetNCells();
1367 Double_t clusE = 0;
1368 UShort_t ids[100] = {0};
1369 Double_t fra[100] = {0};
1370 for (Int_t j = 0; j<nc; ++j) {
1371 Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1372 Double_t cen = cells->GetCellAmplitude(id);
1373 clusE += cen;
1374 if (cen>0) {
1375 ids[nc] = id;
1376 ++nc;
1377 }
1378 }
1379 if (clusE<=0) {
1380 clusters->RemoveAt(i);
1381 continue;
1382 }
1383
1384 for (Int_t j = 0; j<nc; ++j) {
1385 Short_t id = ids[j];
1386 Double_t cen = cells->GetCellAmplitude(id);
1387 fra[j] = cen/clusE;
1388 }
1389 clus->SetE(clusE);
1390 AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1391 if (aodclus) {
1392 aodclus->Clear("");
1393 aodclus->SetNCells(nc);
1394 aodclus->SetCellsAmplitudeFraction(fra);
1395 aodclus->SetCellsAbsId(ids);
1396 }
1397 }
1398 clusters->Compress();
1399}
1400
286b47a5 1401//________________________________________________________________________
1402void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1403{
1404 // Fill histograms related to cell properties.
d595acbb 1405
90d5b88b 1406 AliVCaloCells *cells = fEsdCells;
1407 if (!cells)
1408 cells = fAodCells;
1409
296ea9b4 1410 if (!cells)
1411 return;
90d5b88b 1412
296ea9b4 1413 Int_t cellModCount[12] = {0};
1414 Double_t cellMaxE = 0;
1415 Double_t cellMeanE = 0;
1416 Int_t ncells = cells->GetNumberOfCells();
1417 if (ncells==0)
1418 return;
1419
1420 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1421
1422 for (Int_t i = 0; i<ncells; ++i) {
1423 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1424 Double_t cellE = cells->GetAmplitude(i);
1425 fHCellE->Fill(cellE);
1426 if (cellE>cellMaxE)
1427 cellMaxE = cellE;
1428 cellMeanE += cellE;
1429
1430 Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1431 Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1432 if (!ret) {
1433 AliError(Form("Could not get cell index for %d", absID));
1434 continue;
1435 }
1436 ++cellModCount[iSM];
1437 Int_t iPhi=-1, iEta=-1;
1438 fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1439 fHColuRow[iSM]->Fill(iEta,iPhi,1);
1440 fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1441 fHCellFreqNoCut[iSM]->Fill(absID);
2e4d8148 1442 if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1443 if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
296ea9b4 1444 if (fHCellCheckE && fHCellCheckE[absID])
1445 fHCellCheckE[absID]->Fill(cellE);
2e4d8148 1446 fHCellFreqE[iSM]->Fill(absID, cellE);
296ea9b4 1447 }
1448 fHCellH->Fill(cellMaxE);
1449 cellMeanE /= ncells;
1450 fHCellM->Fill(cellMeanE);
1451 fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1452 for (Int_t i=0; i<nsm; ++i)
1453 fHCellMult[i]->Fill(cellModCount[i]);
286b47a5 1454}
1455
1456//________________________________________________________________________
1457void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1458{
90d5b88b 1459 // Fill histograms related to cluster properties.
fa443410 1460
90d5b88b 1461 TObjArray *clusters = fEsdClusters;
1462 if (!clusters)
1463 clusters = fAodClusters;
296ea9b4 1464 if (!clusters)
1465 return;
90d5b88b 1466
296ea9b4 1467 Double_t vertex[3] = {0};
1468 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1469
1470 Int_t nclus = clusters->GetEntries();
1471 for(Int_t i = 0; i<nclus; ++i) {
1472 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1473 if (!clus)
1474 continue;
1475 if (!clus->IsEMCAL())
1476 continue;
90d5b88b 1477 TLorentzVector clusterVec;
296ea9b4 1478 clus->GetMomentum(clusterVec,vertex);
3a952328 1479 Double_t maxAxis = clus->GetTOF(); //sigma
1480 Double_t clusterEcc = clus->Chi2(); //eccentricity
296ea9b4 1481 fHClustEccentricity->Fill(clusterEcc);
1482 fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1483 fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1484 fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1485 fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1486 fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
f5e0f1e2 1487 fHClustEnergyNCell->Fill(clus->E(),clus->GetNCells());
788ca675 1488 }
1489}
b3ee6797 1490
38727e64 1491//________________________________________________________________________
1492void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
1493{
1494 // Get Mc truth particle information.
1495
cfd7d5b2 1496 if (!fMcParts)
38727e64 1497 return;
1498
807016ea 1499 fMcParts->Clear();
1500
1501 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1502 Double_t etamin = -0.7;
1503 Double_t etamax = +0.7;
1504 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
1505 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
1506
56fd6cb2 1507 if (fAodEv) {
807016ea 1508 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1509 am->LoadBranch(AliAODMCParticle::StdBranchName());
56fd6cb2 1510 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAodEv->FindListObject(AliAODMCParticle::StdBranchName()));
1511 if (!tca)
1512 return;
1513
1514 Int_t nents = tca->GetEntries();
807016ea 1515 for(int it=0; it<nents; ++it) {
56fd6cb2 1516 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(tca->At(it));
807016ea 1517 part->Print();
56fd6cb2 1518
b5ba4932 1519 // pion or eta meson or direct photon
56fd6cb2 1520 if(part->GetPdgCode() == 111) {
1521 } else if(part->GetPdgCode() == 221) {
b5ba4932 1522 } else if(part->GetPdgCode() == 22 ) {
1523 } else
56fd6cb2 1524 continue;
1525
1526 // primary particle
1527 Double_t dR = TMath::Sqrt((part->Xv()*part->Xv())+(part->Yv()*part->Yv()));
1528 if(dR > 1.0)
1529 continue;
56fd6cb2 1530
807016ea 1531 // kinematic cuts
1532 Double_t pt = part->Pt() ;
1533 if (pt<0.5)
1534 continue;
1535 Double_t eta = part->Eta();
1536 if (eta<etamin||eta>etamax)
1537 continue;
1538 Double_t phi = part->Phi();
1539 if (phi<phimin||phi>phimax)
1540 continue;
1541
1542 ProcessDaughters(part, it, tca);
56fd6cb2 1543 }
807016ea 1544 return;
1545 }
38727e64 1546
1547 AliMCEvent *mcEvent = MCEvent();
1548 if (!mcEvent)
1549 return;
1550
1551 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
1552 if (!evtVtx)
1553 return;
1554
1555 mcEvent->PreReadAll();
38727e64 1556
1557 Int_t nTracks = mcEvent->GetNumberOfPrimaries();
1558 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
1559 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1560 if (!mcP)
1561 continue;
1562
b5ba4932 1563 // pion or eta meson or direct photon
807016ea 1564 if(mcP->PdgCode() == 111) {
1565 } else if(mcP->PdgCode() == 221) {
b5ba4932 1566 } else if(mcP->PdgCode() == 22 ) {
38727e64 1567 } else
1568 continue;
1569
1570 // primary particle
807016ea 1571 Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1572 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
38727e64 1573 if(dR > 1.0)
1574 continue;
1575
38727e64 1576 // kinematic cuts
807016ea 1577 Double_t pt = mcP->Pt() ;
1578 if (pt<0.5)
38727e64 1579 continue;
807016ea 1580 Double_t eta = mcP->Eta();
1581 if (eta<etamin||eta>etamax)
38727e64 1582 continue;
807016ea 1583 Double_t phi = mcP->Phi();
1584 if (phi<phimin||phi>phimax)
1585 continue;
1586
1587 ProcessDaughters(mcP, iTrack, mcEvent);
38727e64 1588 }
1589}
1590
788ca675 1591//________________________________________________________________________
1592void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1593{
1594 // Fill ntuple.
1595
1596 if (!fNtuple)
1597 return;
1598
1599 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1600 if (fAodEv) {
1601 fHeader->fRun = fAodEv->GetRunNumber();
1602 fHeader->fOrbit = fAodEv->GetHeader()->GetOrbitNumber();
1603 fHeader->fPeriod = fAodEv->GetHeader()->GetPeriodNumber();
1604 fHeader->fBx = fAodEv->GetHeader()->GetBunchCrossNumber();
1605 fHeader->fL0 = fAodEv->GetHeader()->GetL0TriggerInputs();
1606 fHeader->fL1 = fAodEv->GetHeader()->GetL1TriggerInputs();
1607 fHeader->fL2 = fAodEv->GetHeader()->GetL2TriggerInputs();
1608 fHeader->fTrClassMask = fAodEv->GetHeader()->GetTriggerMask();
1609 fHeader->fTrCluster = fAodEv->GetHeader()->GetTriggerCluster();
1610 fHeader->fOffTriggers = fAodEv->GetHeader()->GetOfflineTrigger();
1611 fHeader->fFiredTriggers = fAodEv->GetHeader()->GetFiredTriggerClasses();
1612 } else {
1613 fHeader->fRun = fEsdEv->GetRunNumber();
1614 fHeader->fOrbit = fEsdEv->GetHeader()->GetOrbitNumber();
1615 fHeader->fPeriod = fEsdEv->GetESDRun()->GetPeriodNumber();
1616 fHeader->fBx = fEsdEv->GetHeader()->GetBunchCrossNumber();
1617 fHeader->fL0 = fEsdEv->GetHeader()->GetL0TriggerInputs();
1618 fHeader->fL1 = fEsdEv->GetHeader()->GetL1TriggerInputs();
1619 fHeader->fL2 = fEsdEv->GetHeader()->GetL2TriggerInputs();
1620 fHeader->fTrClassMask = fEsdEv->GetHeader()->GetTriggerMask();
1621 fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
1622 fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1623 fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
5fe1ca23 1624 Float_t v0CorrR = 0;
1625 fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1626 const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1627 if (mult)
1628 fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1629 fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
c0d2e1f2 1630 AliTriggerAnalysis trAn; /// Trigger Analysis
1631 Bool_t v0B = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0C);
1632 Bool_t v0A = trAn.IsOfflineTriggerFired(fEsdEv, AliTriggerAnalysis::kV0A);
469b2bff 1633 fHeader->fV0And = v0A && v0B;
1634 fHeader->fIsHT = (fHeader->fOffTriggers & AliVEvent::kEMC1) || (fHeader->fOffTriggers & AliVEvent::kEMC7);
1635 am->LoadBranch("SPDPileupVertices");
1636 am->LoadBranch("TrkPileupVertices");
1637 fHeader->fIsPileup = fEsdEv->IsPileupFromSPD(3,0.8);
1638 fHeader->fIsPileup2 = fEsdEv->IsPileupFromSPD(3,0.4);
1639 fHeader->fIsPileup4 = fEsdEv->IsPileupFromSPD(3,0.2);
1640 fHeader->fIsPileup8 = fEsdEv->IsPileupFromSPD(3,0.1);
1641 fHeader->fNSpdVertices = fEsdEv->GetNumberOfPileupVerticesSPD();
1642 fHeader->fNTpcVertices = fEsdEv->GetNumberOfPileupVerticesTracks();
788ca675 1643 }
2ee7bda4 1644
788ca675 1645 AliCentrality *cent = InputEvent()->GetCentrality();
1646 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
1647 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
1648 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
1649 fHeader->fCqual = cent->GetQuality();
b6c599fe 1650
1651 AliEventplane *ep = InputEvent()->GetEventplane();
1652 if (ep) {
1653 if (ep->GetQVector())
1654 fHeader->fPsi = ep->GetQVector()->Phi()/2. ;
f2b8fca6 1655 else
1656 fHeader->fPsi = -1;
b6c599fe 1657 if (ep->GetQsub1()&&ep->GetQsub2())
1658 fHeader->fPsiRes = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1659 else
1660 fHeader->fPsiRes = 0;
1661 }
1662
788ca675 1663 Double_t val = 0;
1664 TString trgclasses(fHeader->fFiredTriggers);
1665 for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1666 const char *name = fTrClassNamesArr->At(j)->GetName();
2ee7bda4 1667 TRegexp regexp(name);
1668 if (trgclasses.Contains(regexp))
788ca675 1669 val += TMath::Power(2,j);
1670 }
1671 fHeader->fTcls = (UInt_t)val;
1672
5fe1ca23 1673 fHeader->fNSelTr = fSelTracks->GetEntries();
1674 fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
f5e0f1e2 1675 fHeader->fNSelPrimTr1 = 0;
1676 fHeader->fNSelPrimTr2 = 0;
1677 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++){
1678 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1679 if(track->Pt()>1)
1680 ++fHeader->fNSelPrimTr1;
1681 if(track->Pt()>2)
1682 ++fHeader->fNSelPrimTr2;
1683 }
5fe1ca23 1684
1685 fHeader->fNCells = 0;
c0d2e1f2 1686 fHeader->fNCells0 = 0;
1687 fHeader->fNCells01 = 0;
1688 fHeader->fNCells03 = 0;
5fe1ca23 1689 fHeader->fNCells1 = 0;
1690 fHeader->fNCells2 = 0;
1691 fHeader->fNCells5 = 0;
1692 fHeader->fMaxCellE = 0;
1693
1694 AliVCaloCells *cells = fEsdCells;
1695 if (!cells)
1696 cells = fAodCells;
1697
1698 if (cells) {
1699 Int_t ncells = cells->GetNumberOfCells();
1700 for(Int_t j=0; j<ncells; ++j) {
1701 Double_t cellen = cells->GetAmplitude(j);
c0d2e1f2 1702 if (cellen>0.045)
1703 ++fHeader->fNCells0;
1704 if (cellen>0.1)
1705 ++fHeader->fNCells01;
1706 if (cellen>0.3)
1707 ++fHeader->fNCells03;
5fe1ca23 1708 if (cellen>1)
1709 ++fHeader->fNCells1;
1710 if (cellen>2)
1711 ++fHeader->fNCells2;
1712 if (cellen>5)
1713 ++fHeader->fNCells5;
1714 if (cellen>fHeader->fMaxCellE)
1715 fHeader->fMaxCellE = cellen;
1716 }
1717 fHeader->fNCells = ncells;
1718 }
1719
1720 fHeader->fNClus = 0;
1721 fHeader->fNClus1 = 0;
1722 fHeader->fNClus2 = 0;
1723 fHeader->fNClus5 = 0;
1724 fHeader->fMaxClusE = 0;
1725
1726 TObjArray *clusters = fEsdClusters;
1727 if (!clusters)
1728 clusters = fAodClusters;
1729
1730 if (clusters) {
1731 Int_t nclus = clusters->GetEntries();
1732 for(Int_t j=0; j<nclus; ++j) {
1733 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
6d4faab0 1734 if (!clus->IsEMCAL())
1735 continue;
5fe1ca23 1736 Double_t clusen = clus->E();
1737 if (clusen>1)
1738 ++fHeader->fNClus1;
1739 if (clusen>2)
1740 ++fHeader->fNClus2;
1741 if (clusen>5)
1742 ++fHeader->fNClus5;
1743 if (clusen>fHeader->fMaxClusE)
1744 fHeader->fMaxClusE = clusen;
1745 }
1746 fHeader->fNClus = nclus;
1747 }
1748
2ee7bda4 1749 fHeader->fMaxTrE = 0;
1750 if (fTriggers) {
1751 Int_t ntrig = fTriggers->GetEntries();
1752 for (Int_t j = 0; j<ntrig; ++j) {
1753 AliStaTrigger *sta = static_cast<AliStaTrigger*>(fTriggers->At(j));
1754 if (!sta)
1755 continue;
1756 if (sta->fE>fHeader->fMaxTrE)
1757 fHeader->fMaxTrE = sta->fE;
1758 }
1759 }
1760
1761 // count cells above 100 MeV on super modules
1762 fHeader->fNcSM0 = GetNCells(0, 0.100);
1763 fHeader->fNcSM1 = GetNCells(1, 0.100);
1764 fHeader->fNcSM2 = GetNCells(2, 0.100);
1765 fHeader->fNcSM3 = GetNCells(3, 0.100);
1766 fHeader->fNcSM4 = GetNCells(4, 0.100);
1767 fHeader->fNcSM5 = GetNCells(5, 0.100);
1768 fHeader->fNcSM6 = GetNCells(6, 0.100);
1769 fHeader->fNcSM7 = GetNCells(7, 0.100);
1770 fHeader->fNcSM8 = GetNCells(8, 0.100);
1771 fHeader->fNcSM9 = GetNCells(9, 0.100);
1772
788ca675 1773 if (fAodEv) {
1774 am->LoadBranch("vertices");
1775 AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1776 FillVertex(fPrimVert, pv);
1777 AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1778 FillVertex(fSpdVert, sv);
1779 } else {
1780 am->LoadBranch("PrimaryVertex.");
1781 const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1782 FillVertex(fPrimVert, pv);
1783 am->LoadBranch("SPDVertex.");
1784 const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1785 FillVertex(fSpdVert, sv);
1786 am->LoadBranch("TPCVertex.");
1787 const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1788 FillVertex(fTpcVert, tv);
fa443410 1789 }
788ca675 1790
788ca675 1791 fNtuple->Fill();
286b47a5 1792}
1793
1794//________________________________________________________________________
1795void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1796{
1797 // Fill histograms related to pions.
286b47a5 1798
296ea9b4 1799 Double_t vertex[3] = {0};
fa443410 1800 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1801
90d5b88b 1802 TObjArray *clusters = fEsdClusters;
1803 if (!clusters)
1804 clusters = fAodClusters;
fa443410 1805
90d5b88b 1806 if (clusters) {
1807 TLorentzVector clusterVec1;
1808 TLorentzVector clusterVec2;
1809 TLorentzVector pionVec;
1810
1811 Int_t nclus = clusters->GetEntries();
fa443410 1812 for (Int_t i = 0; i<nclus; ++i) {
90d5b88b 1813 AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
fa443410 1814 if (!clus1)
1815 continue;
1816 if (!clus1->IsEMCAL())
1817 continue;
3c661da5 1818 if (clus1->E()<fMinE)
6eb6260e 1819 continue;
a49742b5 1820 if (clus1->GetNCells()<fNminCells)
1821 continue;
f224d35b 1822 if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1823 continue;
3c661da5 1824 if (clus1->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1825 continue;
fa443410 1826 clus1->GetMomentum(clusterVec1,vertex);
1827 for (Int_t j = i+1; j<nclus; ++j) {
90d5b88b 1828 AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
fa443410 1829 if (!clus2)
1830 continue;
1831 if (!clus2->IsEMCAL())
1832 continue;
3c661da5 1833 if (clus2->E()<fMinE)
6eb6260e 1834 continue;
a49742b5 1835 if (clus2->GetNCells()<fNminCells)
1836 continue;
f224d35b 1837 if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1838 continue;
3c661da5 1839 if (clus2->Chi2()<fMinEcc) // eccentricity cut
f224d35b 1840 continue;
fa443410 1841 clus2->GetMomentum(clusterVec2,vertex);
1842 pionVec = clusterVec1 + clusterVec2;
1843 Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
6eb6260e 1844 Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
d595acbb 1845 if (pionZgg < fAsymMax) {
1846 fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1847 fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1848 fHPionMggAsym->Fill(pionVec.M(),pionZgg);
6eb6260e 1849 fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
d595acbb 1850 Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1851 fHPionInvMasses[bin]->Fill(pionVec.M());
1852 }
fa443410 1853 }
1854 }
90d5b88b 1855 }
fa443410 1856}
d595acbb 1857
38727e64 1858//________________________________________________________________________
1859void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
1860{
1861 // Fill additional MC information histograms.
1862
cfd7d5b2 1863 if (!fMcParts)
38727e64 1864 return;
1865
1866 // check if aod or esd mc mode and the fill histos
1867}
1868
6eb6260e 1869//________________________________________________________________________
323834f0 1870void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
6eb6260e 1871{
788ca675 1872 // Fill other histograms.
b5ba4932 1873
1874 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); ++iTracks){
1875 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1876 if(!track)
1877 continue;
1878 fHPrimTrackPt->Fill(track->Pt());
1879 fHPrimTrackEta->Fill(track->Eta());
1880 fHPrimTrackPhi->Fill(track->Phi());
1881 }
788ca675 1882}
1883
f5e0f1e2 1884//________________________________________________________________________
1885void AliAnalysisTaskEMCALPi0PbPb::FillTrackHists()
1886{
1887 // Fill track histograms.
1888
1889 if (fSelPrimTracks) {
1890 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++) {
1891 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1892 if(!track)
1893 continue;
1894 fHPrimTrackPt->Fill(track->Pt());
1895 fHPrimTrackEta->Fill(track->Eta());
1896 fHPrimTrackPhi->Fill(track->Phi());
1897 }
1898 }
1899}
1900
788ca675 1901//__________________________________________________________________________________________________
1902void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1903{
1904 // Fill vertex from ESD vertex info.
1905
1906 v->fVx = esdv->GetX();
1907 v->fVy = esdv->GetY();
1908 v->fVz = esdv->GetZ();
1909 v->fVc = esdv->GetNContributors();
1910 v->fDisp = esdv->GetDispersion();
1911 v->fZres = esdv->GetZRes();
1912 v->fChi2 = esdv->GetChi2();
1913 v->fSt = esdv->GetStatus();
1914 v->fIs3D = esdv->IsFromVertexer3D();
1915 v->fIsZ = esdv->IsFromVertexerZ();
1916}
1917
1918//__________________________________________________________________________________________________
1919void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1920{
1921 // Fill vertex from AOD vertex info.
1922
1923 v->fVx = aodv->GetX();
1924 v->fVy = aodv->GetY();
1925 v->fVz = aodv->GetZ();
1926 v->fVc = aodv->GetNContributors();
1927 v->fChi2 = aodv->GetChi2();
6eb6260e 1928}
1929
d595acbb 1930//________________________________________________________________________
296ea9b4 1931Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1932{
1933 // Compute isolation based on cell content.
1934
1935 AliVCaloCells *cells = fEsdCells;
1936 if (!cells)
1937 cells = fAodCells;
1938 if (!cells)
1939 return 0;
1940
1941 Double_t cellIsolation = 0;
1942 Double_t rad2 = radius*radius;
1943 Int_t ncells = cells->GetNumberOfCells();
1944 for (Int_t i = 0; i<ncells; ++i) {
1945 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
0fbe8d4f 1946 Float_t eta=-1, phi=-1;
296ea9b4 1947 fGeom->EtaPhiFromIndex(absID,eta,phi);
0fbe8d4f 1948 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 1949 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1950 if(dist>rad2)
1951 continue;
0fbe8d4f 1952 Double_t cellE = cells->GetAmplitude(i);
2ee7bda4 1953 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
1954 Double_t cellEt = cellE*sin(theta);
1955 cellIsolation += cellEt;
296ea9b4 1956 }
1957 return cellIsolation;
1958}
2ee7bda4 1959
5fc7508c 1960//________________________________________________________________________
1961Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsoNxM(Double_t cEta, Double_t cPhi, Int_t N, Int_t M) const
1962{
1963 // Compute isolation based on cell content, in a NxM rectangle.
1964
1965 AliVCaloCells *cells = fEsdCells;
1966 if (!cells)
1967 cells = fAodCells;
1968 if (!cells)
1969 return 0;
1970
1971 Double_t cellIsolation = 0;
1972 Int_t ncells = cells->GetNumberOfCells();
1973 for (Int_t i = 0; i<ncells; ++i) {
1974 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1975 Float_t eta=-1, phi=-1;
1976 fGeom->EtaPhiFromIndex(absID,eta,phi);
1977 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1978 Double_t etadiff = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1979 if(TMath::Abs(etadiff)/0.014>N)
1980 continue;
1981 if(TMath::Abs(phidiff)/0.014>M)
1982 continue;
1983 Double_t cellE = cells->GetAmplitude(i);
2ee7bda4 1984 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
1985 Double_t cellEt = cellE*sin(theta);
1986 cellIsolation += cellEt;
5fc7508c 1987 }
1988 return cellIsolation;
1989}
1990
296ea9b4 1991//________________________________________________________________________
0fbe8d4f 1992Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
1993{
1994 // Get maximum energy of attached cell.
1995
1996 Double_t ret = 0;
1997 Int_t ncells = cluster->GetNCells();
1998 if (fEsdCells) {
1999 for (Int_t i=0; i<ncells; i++) {
2000 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2001 ret += e;
2002 }
2003 } else {
2004 for (Int_t i=0; i<ncells; i++) {
2005 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2006 ret += e;
2007 }
2008 }
2009 return ret;
2010}
2011
2012//________________________________________________________________________
2013Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
d595acbb 2014{
90d5b88b 2015 // Get maximum energy of attached cell.
2016
788ca675 2017 id = -1;
90d5b88b 2018 Double_t maxe = 0;
90d5b88b 2019 Int_t ncells = cluster->GetNCells();
2020 if (fEsdCells) {
2021 for (Int_t i=0; i<ncells; i++) {
27422847 2022 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
f0897c18 2023 if (e>maxe) {
90d5b88b 2024 maxe = e;
788ca675 2025 id = cluster->GetCellAbsId(i);
f0897c18 2026 }
90d5b88b 2027 }
2028 } else {
2029 for (Int_t i=0; i<ncells; i++) {
27422847 2030 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
90d5b88b 2031 if (e>maxe)
2032 maxe = e;
788ca675 2033 id = cluster->GetCellAbsId(i);
90d5b88b 2034 }
2035 }
6eb6260e 2036 return maxe;
d595acbb 2037}
2038
469b2bff 2039//________________________________________________________________________
1f41ed3d 2040Double_t AliAnalysisTaskEMCALPi0PbPb::GetSecondMaxCellEnergy(AliVCluster *clus, Short_t &id) const
469b2bff 2041{
2042 // Get second maximum cell.
2043
2044 AliVCaloCells *cells = fEsdCells;
2045 if (!cells)
2046 cells = fAodCells;
2047 if (!cells)
2048 return -1;
2049
2050 Double_t secondEmax=0, firstEmax=0;
2051 Double_t cellen;
2052 for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2053 Int_t absId = clus->GetCellAbsId(iCell);
2054 cellen = cells->GetCellAmplitude(absId);
2055 if(cellen > firstEmax)
2056 firstEmax = cellen;
2057 }
2058 for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
2059 Int_t absId = clus->GetCellAbsId(iCell);
2060 cellen = cells->GetCellAmplitude(absId);
1f41ed3d 2061 if(cellen < firstEmax && cellen > secondEmax) {
469b2bff 2062 secondEmax = cellen;
1f41ed3d 2063 id = absId;
2064 }
469b2bff 2065 }
2066 return secondEmax;
2067}
2068
d595acbb 2069//________________________________________________________________________
0fbe8d4f 2070void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
d595acbb 2071{
90d5b88b 2072 // Calculate the (E) weighted variance along the longer (eigen) axis.
2073
6bf90832 2074 sigmaMax = 0; // cluster variance along its longer axis
2075 sigmaMin = 0; // cluster variance along its shorter axis
2076 Double_t Ec = c->E(); // cluster energy
296ea9b4 2077 if(Ec<=0)
2078 return;
6bf90832 2079 Double_t Xc = 0 ; // cluster first moment along X
2080 Double_t Yc = 0 ; // cluster first moment along Y
2081 Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2)
2082 Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2)
2083 Double_t Syy = 0 ; // cluster covariance^2
90d5b88b 2084
2085 AliVCaloCells *cells = fEsdCells;
2086 if (!cells)
2087 cells = fAodCells;
2088
6eb6260e 2089 if (!cells)
2090 return;
2091
6bf90832 2092 Int_t ncells = c->GetNCells();
6eb6260e 2093 if (ncells==1)
2094 return;
2095
6eb6260e 2096 for(Int_t j=0; j<ncells; ++j) {
6bf90832 2097 Int_t id = TMath::Abs(c->GetCellAbsId(j));
6eb6260e 2098 Double_t cellen = cells->GetCellAmplitude(id);
3a952328 2099 TVector3 pos;
2100 fGeom->GetGlobal(id,pos);
6eb6260e 2101 Xc += cellen*pos.X();
2102 Yc += cellen*pos.Y();
2103 Sxx += cellen*pos.X()*pos.X();
2104 Syy += cellen*pos.Y()*pos.Y();
2105 Sxy += cellen*pos.X()*pos.Y();
2106 }
2107 Xc /= Ec;
2108 Yc /= Ec;
2109 Sxx /= Ec;
2110 Syy /= Ec;
2111 Sxy /= Ec;
2112 Sxx -= Xc*Xc;
2113 Syy -= Yc*Yc;
2114 Sxy -= Xc*Yc;
f0897c18 2115 Sxx = TMath::Abs(Sxx);
2116 Syy = TMath::Abs(Syy);
296ea9b4 2117 sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2118 sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
2119 sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2120 sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
d595acbb 2121}
2ee7bda4 2122
5fc7508c 2123//________________________________________________________________________
2124void AliAnalysisTaskEMCALPi0PbPb::GetSigmaEtaEta(const AliVCluster *c, Double_t& sEtaEta, Double_t &sPhiPhi) const
2125{
2126 // Calculate the (E) weighted variance along the pseudorapidity.
2127
2128 sEtaEta = 0;
2129 sPhiPhi = 0;
2130
2131 Double_t Ec = c->E(); // cluster energy
2132 if(Ec<=0)
2133 return;
2134
2135 const Int_t ncells = c->GetNCells();
2136
2137 Double_t EtaC = 0; // cluster first moment along eta
2138 Double_t PhiC = 0; // cluster first moment along phi
2139 Double_t Setaeta = 0; // cluster second central moment along eta
2140 Double_t Sphiphi = 0; // cluster second central moment along phi
2141 Double_t w[ncells]; // weight max(0,4.5*log(E_i/Ec))
2142 Double_t sumw = 0;
2143 Int_t id[ncells];
2144
2145 AliVCaloCells *cells = fEsdCells;
2146 if (!cells)
2147 cells = fAodCells;
2148
2149 if (!cells)
2150 return;
2151
2152 if (ncells==1)
2153 return;
2154
2155 for(Int_t j=0; j<ncells; ++j) {
2156 id[j] = TMath::Abs(c->GetCellAbsId(j));
2157 Double_t cellen = cells->GetCellAmplitude(id[j]);
2158 w[j] = TMath::Max(0., 4.5+TMath::Log(cellen/Ec));
2159 TVector3 pos;
2160 fGeom->GetGlobal(id[j],pos);
2161 EtaC += w[j]*pos.Eta();
2162 PhiC += w[j]*pos.Phi();
2163 sumw += w[j];
2164 }
2165 EtaC /= sumw;
2166 PhiC /= sumw;
2167
2168 for(Int_t j=0; j<ncells; ++j) {
2169 TVector3 pos;
2170 fGeom->GetGlobal(id[j],pos);
2171 Setaeta = w[j]*(pos.Eta() - EtaC)*(pos.Eta() - EtaC);
2172 Sphiphi = w[j]*(pos.Phi() - PhiC)*(pos.Phi() - PhiC);
2173 }
2174 Setaeta /= sumw;
2175 sEtaEta = TMath::Sqrt(Setaeta);
2176 Sphiphi /= sumw;
2177 sPhiPhi = TMath::Sqrt(Sphiphi);
2178}
2179
6bf90832 2180//________________________________________________________________________
0fbe8d4f 2181Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
6bf90832 2182{
2183 // Calculate number of attached cells above emin.
2184
6bf90832 2185 AliVCaloCells *cells = fEsdCells;
2186 if (!cells)
2187 cells = fAodCells;
6bf90832 2188 if (!cells)
296ea9b4 2189 return 0;
6bf90832 2190
296ea9b4 2191 Int_t n = 0;
6bf90832 2192 Int_t ncells = c->GetNCells();
2193 for(Int_t j=0; j<ncells; ++j) {
2194 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2195 Double_t cellen = cells->GetCellAmplitude(id);
2196 if (cellen>=emin)
2197 ++n;
2198 }
2199 return n;
2200}
2201
2ee7bda4 2202//________________________________________________________________________
2203Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(Int_t sm, Double_t emin) const
2204{
2205 // Calculate number of cells per SM above emin.
2206
2207 AliVCaloCells *cells = fEsdCells;
2208 if (!cells)
2209 cells = fAodCells;
2210 if (!cells)
2211 return 0;
2212
2213 Int_t n = 0;
2214 Int_t ncells = cells->GetNumberOfCells();
2215 for(Int_t j=0; j<ncells; ++j) {
2216 Int_t id = TMath::Abs(cells->GetCellNumber(j));
2217 Double_t cellen = cells->GetCellAmplitude(id);
2218 if (cellen<emin)
2219 continue;
2220 Int_t fsm = fGeom->GetSuperModuleNumber(id);
2221 if (fsm != sm)
2222 continue;
2223 ++n;
2224 }
2225 return n;
2226}
2227
296ea9b4 2228//________________________________________________________________________
b6c599fe 2229Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
296ea9b4 2230{
2231 // Compute isolation based on tracks.
2232
2233 Double_t trkIsolation = 0;
2234 Double_t rad2 = radius*radius;
b6c599fe 2235 Int_t ntrks = fSelPrimTracks->GetEntries();
296ea9b4 2236 for(Int_t j = 0; j<ntrks; ++j) {
1f41ed3d 2237 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
296ea9b4 2238 if (!track)
2239 continue;
b6c599fe 2240 if (track->Pt()<pt)
2241 continue;
296ea9b4 2242 Float_t eta = track->Eta();
2243 Float_t phi = track->Phi();
0fbe8d4f 2244 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
296ea9b4 2245 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
2246 if(dist>rad2)
2247 continue;
2248 trkIsolation += track->Pt();
2249 }
2250 return trkIsolation;
2251}
2ee7bda4 2252
5fc7508c 2253//________________________________________________________________________
2254Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsoStrip(Double_t cEta, Double_t cPhi, Double_t dEta, Double_t dPhi, Double_t pt) const
2255{
2256 // Compute isolation based on tracks.
2257
2258 Double_t trkIsolation = 0;
2259 Int_t ntrks = fSelPrimTracks->GetEntries();
2260 for(Int_t j = 0; j<ntrks; ++j) {
1f41ed3d 2261 AliVTrack *track = static_cast<AliVTrack*>(fSelPrimTracks->At(j));
5fc7508c 2262 if (!track)
2263 continue;
2264 if (track->Pt()<pt)
2265 continue;
2266 Float_t eta = track->Eta();
2267 Float_t phi = track->Phi();
2268 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2269 Double_t etadiff = (eta-cEta);
2270 if(TMath::Abs(etadiff)>dEta || TMath::Abs(phidiff)>dPhi)
2271 continue;
2272 trkIsolation += track->Pt();
2273 }
2274 return trkIsolation;
2275}
2276
0fbe8d4f 2277//________________________________________________________________________
2278Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
2279{
2280 // Returns if cluster shared across super modules.
2281
1f41ed3d 2282 if (!c)
0fbe8d4f 2283 return 0;
2284
2285 Int_t n = -1;
2286 Int_t ncells = c->GetNCells();
2287 for(Int_t j=0; j<ncells; ++j) {
2288 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2289 Int_t got = id / (24*48);
2290 if (n==-1) {
2291 n = got;
2292 continue;
2293 }
2294 if (got!=n)
2295 return 1;
2296 }
2297 return 0;
2298}
3a952328 2299
2ee7bda4 2300//________________________________________________________________________
2301Bool_t AliAnalysisTaskEMCALPi0PbPb::IsIdPartOfCluster(const AliVCluster *c, Short_t id) const
2302{
2303 // Returns if id is part of cluster.
2304
2305 AliVCaloCells *cells = fEsdCells;
2306 if (!cells)
2307 cells = fAodCells;
2308 if (!cells)
2309 return 0;
2310
2311 Int_t ncells = c->GetNCells();
2312 for(Int_t j=0; j<ncells; ++j) {
2313 Int_t cid = TMath::Abs(c->GetCellAbsId(j));
2314 if (cid == id)
2315 return 1;
2316 }
2317 return 0;
2318}
2319
38727e64 2320//________________________________________________________________________
807016ea 2321void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliVParticle *p, const TObjArray *arr, Int_t level) const
38727e64 2322{
56fd6cb2 2323 // Print recursively daughter information.
2324
2325 if (!p || !arr)
2326 return;
2327
807016ea 2328 const AliAODMCParticle *amc = dynamic_cast<const AliAODMCParticle*>(p);
2329 if (!amc)
2330 return;
2331 for (Int_t i=0; i<level; ++i) printf(" ");
2332 amc->Print();
2333
2334 Int_t n = amc->GetNDaughters();
2335 for (Int_t i=0; i<n; ++i) {
2336 Int_t d = amc->GetDaughter(i);
2337 const AliVParticle *dmc = static_cast<const AliVParticle*>(arr->At(d));
2338 PrintDaughters(dmc,arr,level+1);
2339 }
2340}
2341
2342//________________________________________________________________________
2343void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliMCParticle *p, const AliMCEvent *arr, Int_t level) const
2344{
2345 // Print recursively daughter information.
2346
2347 if (!p || !arr)
2348 return;
2349
2350 for (Int_t i=0; i<level; ++i) printf(" ");
2351 Int_t d1 = p->GetFirstDaughter();
2352 Int_t d2 = p->GetLastDaughter();
2353 printf("pid=%d: %.2f %.2f %.2f (%.2f %.2f %.2f); nd=%d,%d\n",
2354 p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2);
2355 if (d1<0)
2356 return;
2357 if (d2<0)
2358 d2=d1;
2359 for (Int_t i=d1;i<=d2;++i) {
2360 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(arr->GetTrack(i));
6d4faab0 2361 PrintDaughters(dmc,arr,level+1);
56fd6cb2 2362 }
38727e64 2363}
2364
2365//________________________________________________________________________
2366void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
2367{
56fd6cb2 2368 // Print track ref array.
2369
2370 if (!p)
2371 return;
2372
2373 Int_t n = p->GetNumberOfTrackReferences();
2374 for (Int_t i=0; i<n; ++i) {
2375 AliTrackReference *ref = p->GetTrackReference(i);
bc8a45bd 2376 if (!ref)
2377 continue;
56fd6cb2 2378 ref->SetUserId(ref->DetectorId());
2379 ref->Print();
2380 }
38727e64 2381}
2382
807016ea 2383//________________________________________________________________________
2384void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliVParticle *p, Int_t index, const TObjArray *arr)
2385{
2386 // Process and create daughters.
2387
2388 if (!p || !arr)
2389 return;
2390
2391 AliAODMCParticle *amc = dynamic_cast<AliAODMCParticle*>(p);
2392 if (!amc)
2393 return;
2394
2395 //amc->Print();
2396
2397 Int_t nparts = arr->GetEntries();
2398 Int_t nents = fMcParts->GetEntries();
2399
2400 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2401 newp->fPt = amc->Pt();
2402 newp->fEta = amc->Eta();
2403 newp->fPhi = amc->Phi();
2404 if (amc->Xv() != 0 || amc->Yv() != 0 || amc->Zv() != 0) {
2405 TVector3 vec(amc->Xv(),amc->Yv(),amc->Zv());
2406 newp->fVR = vec.Perp();
2407 newp->fVEta = vec.Eta();
2408 newp->fVPhi = vec.Phi();
2409 }
2410 newp->fPid = amc->PdgCode();
8c56d760 2411 newp->fLab = nents;
807016ea 2412 Int_t moi = amc->GetMother();
2413 if (moi>=0&&moi<nparts) {
2414 const AliAODMCParticle *mmc = static_cast<const AliAODMCParticle*>(arr->At(moi));
2415 moi = mmc->GetUniqueID();
2416 }
2417 newp->fMo = moi;
2418 p->SetUniqueID(nents);
8c56d760 2419
2420 // TODO: Determine which detector was hit
2421 //newp->fDet = ???
2422
807016ea 2423 Int_t n = amc->GetNDaughters();
2424 for (Int_t i=0; i<n; ++i) {
2425 Int_t d = amc->GetDaughter(i);
2426 if (d<=index || d>=nparts)
2427 continue;
2428 AliVParticle *dmc = static_cast<AliVParticle*>(arr->At(d));
2429 ProcessDaughters(dmc,d,arr);
2430 }
2431}
2432
2433//________________________________________________________________________
2434void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliMCParticle *p, Int_t index, const AliMCEvent *arr)
2435{
2436 // Process and create daughters.
2437
2438 if (!p || !arr)
2439 return;
2440
2441 Int_t d1 = p->GetFirstDaughter();
2442 Int_t d2 = p->GetLastDaughter();
2443 if (0) {
2444 printf("%d pid=%d: %.3f %.3f %.3f (%.2f %.2f %.2f); nd=%d,%d, mo=%d\n",
2445 index,p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2, p->GetMother());
bc8a45bd 2446 PrintTrackRefs(p);
807016ea 2447 }
2448 Int_t nents = fMcParts->GetEntries();
2449
2450 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2451 newp->fPt = p->Pt();
2452 newp->fEta = p->Eta();
2453 newp->fPhi = p->Phi();
2454 if (p->Xv() != 0 || p->Yv() != 0 || p->Zv() != 0) {
2455 TVector3 vec(p->Xv(),p->Yv(),p->Zv());
2456 newp->fVR = vec.Perp();
2457 newp->fVEta = vec.Eta();
2458 newp->fVPhi = vec.Phi();
2459 }
2460 newp->fPid = p->PdgCode();
8c56d760 2461 newp->fLab = nents;
807016ea 2462 Int_t moi = p->GetMother();
2463 if (moi>=0) {
2464 const AliMCParticle *mmc = static_cast<const AliMCParticle *>(arr->GetTrack(moi));
2465 moi = mmc->GetUniqueID();
2466 }
2467 newp->fMo = moi;
2468 p->SetUniqueID(nents);
2469
2470 Int_t nref = p->GetNumberOfTrackReferences();
2471 if (nref>0) {
2472 AliTrackReference *ref = p->GetTrackReference(nref-1);
2473 if (ref) {
2474 newp->fDet = ref->DetectorId();
2475 }
2476 }
2477
2478 if (d1<0)
2479 return;
2480 if (d2<0)
2481 d2=d1;
2482 for (Int_t i=d1;i<=d2;++i) {
2483 AliMCParticle *dmc = static_cast<AliMCParticle *>(arr->GetTrack(i));
2484 if (dmc->P()<0.01)
2485 continue;
2486 ProcessDaughters(dmc,i,arr);
2487 }
2488}
6dbf6b27 2489
2490//__________________________________________________________________________________________________
2491void AliStaCluster::GetMom(TLorentzVector& p, Double_t *vertex)
2492{
2493 // Calculate momentum.
2494
2495 TVector3 pos;
2496 pos.SetPtEtaPhi(fR,fEta,fPhi);
2497
2498 if(vertex){ //calculate direction relative to vertex
2499 pos -= vertex;
2500 }
2501
2502 Double_t r = pos.Mag();
2503 p.SetPxPyPzE(fE*pos.x()/r, fE*pos.y()/r, fE*pos.z()/r, fE);
2504}
2505
2506//__________________________________________________________________________________________________
2507void AliStaCluster::GetMom(TLorentzVector& p, AliStaVertex *vertex)
2508{
2509 // Calculate momentum.
2510
2511 Double_t v[3] = {0,0,0};
2512 if (vertex) {
2513 v[0] = vertex->fVx;
2514 v[1] = vertex->fVy;
2515 v[2] = vertex->fVz;
2516 }
2517 GetMom(p, v);
2518}