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