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