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