3 #include "AliAnalysisTaskEMCALPi0PbPb.h"
6 #include <TClonesArray.h>
8 #include <TGeoGlobalMagField.h>
12 #include <TLorentzVector.h>
17 #include "AliAODEvent.h"
18 #include "AliAODMCParticle.h"
19 #include "AliAODVertex.h"
20 #include "AliAnalysisManager.h"
21 #include "AliAnalysisTaskEMCALClusterizeFast.h"
22 #include "AliCDBManager.h"
23 #include "AliCentrality.h"
24 #include "AliEMCALDigit.h"
25 #include "AliEMCALGeometry.h"
26 #include "AliEMCALRecPoint.h"
27 #include "AliEMCALRecoUtils.h"
28 #include "AliESDCaloTrigger.h"
29 #include "AliESDEvent.h"
30 #include "AliESDUtils.h"
31 #include "AliESDVertex.h"
32 #include "AliESDtrack.h"
33 #include "AliESDtrackCuts.h"
34 #include "AliEventplane.h"
35 #include "AliGeomManager.h"
36 #include "AliInputEventHandler.h"
38 #include "AliMCEvent.h"
39 #include "AliMCParticle.h"
41 #include "AliMultiplicity.h"
43 #include "AliTrackerBase.h"
45 ClassImp(AliAnalysisTaskEMCALPi0PbPb)
47 //________________________________________________________________________
48 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb()
49 : AliAnalysisTaskSE(),
64 fGeoName("EMCAL_FIRSTYEARV1"),
110 fHTclsBeforeCuts(0x0),
111 fHTclsAfterCuts(0x0),
119 fHCellFreqNoCut(0x0),
120 fHCellFreqCut100M(0x0),
121 fHCellFreqCut300M(0x0),
124 fHClustEccentricity(0),
126 fHClustEnergyPt(0x0),
127 fHClustEnergySigma(0x0),
128 fHClustSigmaSigma(0x0),
129 fHClustNCellEnergyRatio(0x0),
130 fHClustEnergyNCell(0x0),
145 //________________________________________________________________________
146 AliAnalysisTaskEMCALPi0PbPb::AliAnalysisTaskEMCALPi0PbPb(const char *name)
147 : AliAnalysisTaskSE(name),
162 fGeoName("EMCAL_FIRSTYEARV1"),
208 fHTclsBeforeCuts(0x0),
209 fHTclsAfterCuts(0x0),
217 fHCellFreqNoCut(0x0),
218 fHCellFreqCut100M(0x0),
219 fHCellFreqCut300M(0x0),
222 fHClustEccentricity(0),
224 fHClustEnergyPt(0x0),
225 fHClustEnergySigma(0x0),
226 fHClustSigmaSigma(0x0),
227 fHClustNCellEnergyRatio(0x0),
228 fHClustEnergyNCell(0x0),
242 DefineOutput(1, TList::Class());
243 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.,SPDVertex.,TPCVertex.,EMCALCells.,Tracks EMCALTrigger."
244 "AOD:header,vertices,emcalCells,tracks";
247 //________________________________________________________________________
248 AliAnalysisTaskEMCALPi0PbPb::~AliAnalysisTaskEMCALPi0PbPb()
252 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
253 delete fOutput; fOutput = 0;
255 delete fPtRanges; fPtRanges = 0;
256 delete fGeom; fGeom = 0;
257 delete fReco; fReco = 0;
258 delete fTrClassNamesArr;
260 delete fSelPrimTracks;
261 delete [] fAmpInTrigger;
263 delete [] fHColuRowE;
264 delete [] fHCellMult;
265 delete [] fHCellFreqNoCut;
266 delete [] fHCellFreqCut100M;
267 delete [] fHCellFreqCut300M;
270 //________________________________________________________________________
271 void AliAnalysisTaskEMCALPi0PbPb::UserCreateOutputObjects()
273 // Create user objects here.
275 cout << "AliAnalysisTaskEMCALPi0PbPb: Input settings" << endl;
276 cout << " fCentVar: " << fCentVar << endl;
277 cout << " fCentFrom: " << fCentFrom << endl;
278 cout << " fCentTo: " << fCentTo << endl;
279 cout << " fVtxZMin: " << fVtxZMin << endl;
280 cout << " fVtxZMax: " << fVtxZMax << endl;
281 cout << " fUseQualFlag: " << fUseQualFlag << endl;
282 cout << " fClusName: \"" << fClusName << "\"" << endl;
283 cout << " fDoNtuple: " << fDoNtuple << endl;
284 cout << " fDoAfterburner: " << fDoAfterburner << endl;
285 cout << " fAsymMax: " << fAsymMax << endl;
286 cout << " fNminCells: " << fNminCells << endl;
287 cout << " fMinE: " << fMinE << endl;
288 cout << " fMinErat: " << fMinErat << endl;
289 cout << " fMinEcc: " << fMinEcc << endl;
290 cout << " fGeoName: \"" << fGeoName << "\"" << endl;
291 cout << " fMinNClusPerTr: " << fMinNClusPerTr << endl;
292 cout << " fIsoDist: " << fIsoDist << endl;
293 cout << " fTrClassNames: \"" << fTrClassNames << "\"" << endl;
294 cout << " fTrCuts: " << fTrCuts << endl;
295 cout << " fPrimTrCuts: " << fPrimTrCuts << endl;
296 cout << " fDoTrMatGeom: " << fDoTrMatGeom << endl;
297 cout << " fTrainMode: " << fTrainMode << endl;
298 cout << " fMarkCells: " << fMarkCells << endl;
299 cout << " fMinL0Time: " << fMinL0Time << endl;
300 cout << " fMaxL0Time: " << fMaxL0Time << endl;
301 cout << " fMcMode: " << fMcMode << endl;
302 cout << " fEmbedMode: " << fEmbedMode << endl;
303 cout << " fGeom: " << fGeom << endl;
304 cout << " fReco: " << fReco << endl;
305 cout << " fDoPSel: " << fDoPSel << endl;
308 fGeom = AliEMCALGeometry::GetInstance(fGeoName);
310 if (fGeom->GetMatrixForSuperModule(0))
311 fIsGeoMatsSet = kTRUE;
314 fReco = new AliEMCALRecoUtils();
315 fTrClassNamesArr = fTrClassNames.Tokenize(" ");
316 fOutput = new TList();
317 fOutput->SetOwner(1);
318 fSelTracks = new TObjArray;
319 fSelPrimTracks = new TObjArray;
322 TFile *f = OpenFile(1);
323 TDirectory::TContext context(f);
325 f->SetCompressionLevel(2);
326 fNtuple = new TTree(Form("tree%.0fto%.0f",fCentFrom,fCentTo), "StandaloneTree");
327 fNtuple->SetDirectory(f);
329 fNtuple->SetAutoFlush(-2*1024*1024);
330 fNtuple->SetAutoSave(0);
332 fNtuple->SetAutoFlush(-32*1024*1024);
333 fNtuple->SetAutoSave(0);
336 fHeader = new AliStaHeader;
337 fNtuple->Branch("header", &fHeader, 16*1024, 99);
338 fPrimVert = new AliStaVertex;
339 fNtuple->Branch("primv", &fPrimVert, 16*1024, 99);
340 fSpdVert = new AliStaVertex;
341 fNtuple->Branch("spdv", &fSpdVert, 16*1024, 99);
342 fTpcVert = new AliStaVertex;
343 fNtuple->Branch("tpcv", &fTpcVert, 16*1024, 99);
344 if (TClass::GetClass("AliStaCluster"))
345 TClass::GetClass("AliStaCluster")->IgnoreTObjectStreamer();
346 fClusters = new TClonesArray("AliStaCluster");
347 fNtuple->Branch("clusters", &fClusters, 8*16*1024, 99);
348 if (TClass::GetClass("AliStaTrigger"))
349 TClass::GetClass("AliStaTrigger")->IgnoreTObjectStreamer();
350 fTriggers = new TClonesArray("AliStaTrigger");
351 fNtuple->Branch("l0prim", &fTriggers, 16*1024, 99);
352 if (fMcMode||fEmbedMode) {
353 if (TClass::GetClass("AliStaPart"))
354 TClass::GetClass("AliStaPart")->IgnoreTObjectStreamer();
355 fMcParts = new TClonesArray("AliStaPart");
356 fNtuple->Branch("mcparts", &fMcParts, 8*16*1024, 99);
361 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
362 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
363 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
366 TH1::SetDefaultSumw2(kTRUE);
367 TH2::SetDefaultSumw2(kTRUE);
368 fHCuts = new TH1F("hEventCuts","",5,0.5,5.5);
369 fHCuts->GetXaxis()->SetBinLabel(1,"All");
370 fHCuts->GetXaxis()->SetBinLabel(2,"PS");
371 fHCuts->GetXaxis()->SetBinLabel(3,Form("%s: %.0f-%.0f",fCentVar.Data(),fCentFrom,fCentTo));
372 fHCuts->GetXaxis()->SetBinLabel(4,"QFlag");
373 fHCuts->GetXaxis()->SetBinLabel(5,Form("zvtx: %.0f-%.0f",fVtxZMin,fVtxZMax));
374 fOutput->Add(fHCuts);
375 fHVertexZ = new TH1F("hVertexZBeforeCut","",100,-25,25);
376 fHVertexZ->SetXTitle("z [cm]");
377 fOutput->Add(fHVertexZ);
378 fHVertexZ2 = new TH1F("hVertexZAfterCut","",100,-25,25);
379 fHVertexZ2->SetXTitle("z [cm]");
380 fOutput->Add(fHVertexZ2);
381 fHCent = new TH1F("hCentBeforeCut","",102,-1,101);
382 fHCent->SetXTitle(fCentVar.Data());
383 fOutput->Add(fHCent);
384 fHCentQual = new TH1F("hCentAfterCut","",102,-1,101);
385 fHCentQual->SetXTitle(fCentVar.Data());
386 fOutput->Add(fHCentQual);
387 fHTclsBeforeCuts = new TH1F("hTclsBeforeCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
388 fHTclsAfterCuts = new TH1F("hTclsAfterCuts","",fTrClassNamesArr->GetEntries(),0.5,0.5+fTrClassNamesArr->GetEntries());
389 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
390 const char *name = fTrClassNamesArr->At(i)->GetName();
391 fHTclsBeforeCuts->GetXaxis()->SetBinLabel(1+i,name);
392 fHTclsAfterCuts->GetXaxis()->SetBinLabel(1+i,name);
394 fOutput->Add(fHTclsBeforeCuts);
395 fOutput->Add(fHTclsAfterCuts);
397 // histograms for cells
398 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
399 fHColuRow = new TH2*[nsm];
400 fHColuRowE = new TH2*[nsm];
401 fHCellMult = new TH1*[nsm];
402 for (Int_t i = 0; i<nsm; ++i) {
403 fHColuRow[i] = new TH2F(Form("hColRow_Mod%d",i),"",48,0,48,24,0.,24);
404 fHColuRow[i]->SetTitle(Form("Module %d: Occupancy", i));
405 fHColuRow[i]->SetXTitle("col (i#eta)");
406 fHColuRow[i]->SetYTitle("row (i#phi)");
407 fHColuRowE[i] = new TH2F(Form("hColRowE_Mod%d", i),"",48,0,48,24,0,24);
408 fHColuRowE[i]->SetTitle(Form("Module %d: Cell energy",i));
409 fHColuRowE[i]->SetXTitle("col (i#eta)");
410 fHColuRowE[i]->SetYTitle("row (i#phi)");
411 fHCellMult[i] = new TH1F(Form("hCellMult_Mod%d",i),"",1000,0,1000);
412 fHCellMult[i]->SetTitle(Form("Module %d: Cell multiplicity",i));
413 fHCellMult[i]->SetXTitle("# of cells");
414 fOutput->Add(fHColuRow[i]);
415 fOutput->Add(fHColuRowE[i]);
416 fOutput->Add(fHCellMult[i]);
418 fHCellE = new TH1F("hCellE","",250,0.,25.);
419 fHCellE->SetXTitle("E_{cell} [GeV]");
420 fOutput->Add(fHCellE);
421 fHCellH = new TH1F ("hCellHighestE","",250,0.,25.);
422 fHCellH->SetXTitle("E^{max}_{cell} [GeV]");
423 fOutput->Add(fHCellH);
424 fHCellM = new TH1F ("hCellMeanEperHitCell","",250,0.,2.5);
425 fHCellM->SetXTitle("#LT E_{cell}#GT [GeV]");
426 fOutput->Add(fHCellM);
427 fHCellM2 = new TH1F ("hCellMeanEperAllCells","",250,0.,1);
428 fHCellM2->SetXTitle("1/N_{cells} #Sigma E_{cell} [GeV]");
429 fOutput->Add(fHCellM2);
431 fHCellFreqNoCut = new TH1*[nsm];
432 fHCellFreqCut100M = new TH1*[nsm];
433 fHCellFreqCut300M = new TH1*[nsm];
434 fHCellFreqE = new TH1*[nsm];
435 for (Int_t i = 0; i<nsm; ++i){
436 Double_t lbin = i*24*48-0.5;
437 Double_t hbin = lbin+24*48;
438 fHCellFreqNoCut[i] = new TH1F(Form("hCellFreqNoCut_SM%d",i),
439 Form("Frequency SM%d (no cut);id;#",i), 1152, lbin, hbin);
440 fHCellFreqCut100M[i] = new TH1F(Form("hCellFreqCut100M_SM%d",i),
441 Form("Frequency SM%d (>0.1GeV);id;#",i), 1152, lbin, hbin);
442 fHCellFreqCut300M[i] = new TH1F(Form("hCellFreqCut300M_SM%d",i),
443 Form("Frequency SM%d (>0.3GeV);id;#",i), 1152, lbin, hbin);
444 fHCellFreqE[i] = new TH1F(Form("hCellFreqE_SM%d",i),
445 Form("Frequency SM%d (E weighted);id;#",i), 1152, lbin, hbin);
446 fOutput->Add(fHCellFreqNoCut[i]);
447 fOutput->Add(fHCellFreqCut100M[i]);
448 fOutput->Add(fHCellFreqCut300M[i]);
449 fOutput->Add(fHCellFreqE[i]);
451 if (!fMarkCells.IsNull()) {
452 fHCellCheckE = new TH1*[24*48*nsm];
453 memset(fHCellCheckE,0,24*48*nsm*sizeof(TH1*));
454 TObjArray *cells = fMarkCells.Tokenize(" ");
455 Int_t n = cells->GetEntries();
456 Int_t *tcs = new Int_t[n];
457 for (Int_t i=0;i<n;++i) {
458 TString name(cells->At(i)->GetName());
461 for (Int_t i = 0; i<n; ++i) {
464 fHCellCheckE[i] = new TH1F(Form("hCellE_id%d",c), Form("Cell %d;E [GeV/c];#",c), 1000, 0, 10);
465 fOutput->Add(fHCellCheckE[i]);
472 // histograms for clusters
474 fHClustEccentricity = new TH1F("hClustEccentricity","",100,-0.1,1.1);
475 fHClustEccentricity->SetXTitle("#epsilon_{C}");
476 fOutput->Add(fHClustEccentricity);
477 fHClustEtaPhi = new TH2F("hClustEtaPhi","",500,-0.8,0.8,100*nsm,phimin,phimax);
478 fHClustEtaPhi->SetXTitle("#eta");
479 fHClustEtaPhi->SetYTitle("#varphi");
480 fOutput->Add(fHClustEtaPhi);
481 fHClustEnergyPt = new TH2F("hClustEnergyPt","",250,0,50,250,0,50);
482 fHClustEnergyPt->SetXTitle("E [GeV]");
483 fHClustEnergyPt->SetYTitle("p_{T} [GeV/c]");
484 fOutput->Add(fHClustEnergyPt);
485 fHClustEnergySigma = new TH2F("hClustEnergySigma","",100,0,100,500,0,50);
486 fHClustEnergySigma->SetXTitle("E_{C} * #sigma_{max} [GeV*cm]");
487 fHClustEnergySigma->SetYTitle("E_{C} [GeV]");
488 fOutput->Add(fHClustEnergySigma);
489 fHClustSigmaSigma = new TH2F("hClustSigmaSigma","",500,0,50,500,0,50);
490 fHClustSigmaSigma->SetXTitle("#lambda_{0} [cm]");
491 fHClustSigmaSigma->SetYTitle("#sigma_{max} [cm]");
492 fOutput->Add(fHClustSigmaSigma);
493 fHClustNCellEnergyRatio = new TH2F("hClustNCellEnergyRatio","",27,-0.5,26.5,101,-0.05,1.05);
494 fHClustNCellEnergyRatio->SetXTitle("N_{cells}");
495 fHClustNCellEnergyRatio->SetYTitle("E^{max}_{cell}/E_{clus}");
496 fOutput->Add(fHClustNCellEnergyRatio);
497 fHClustEnergyNCell = new TH2F("hClustEnergyNCell","",200,0,100,50,0,50);
498 fHClustEnergyNCell->SetXTitle("E_{clus}");
499 fHClustEnergyNCell->SetYTitle("N_{cells}");
500 fOutput->Add(fHClustEnergyNCell);
503 // histograms for primary tracks
504 fHPrimTrackPt = new TH1F("hPrimTrackPt",";p_{T} [GeV/c]",500,0,50);
505 fOutput->Add(fHPrimTrackPt);
506 fHPrimTrackEta = new TH1F("hPrimTrackEta",";#eta",40,-2,2);
507 fOutput->Add(fHPrimTrackEta);
508 fHPrimTrackPhi = new TH1F("hPrimTrackPhi",";#varPhi [rad]",63,0,6.3);
509 fOutput->Add(fHPrimTrackPhi);
510 // histograms for track matching
511 fHMatchDr = new TH1F("hMatchDrDist",";dR [cm]",500,0,200);
512 fOutput->Add(fHMatchDr);
513 fHMatchDz = new TH1F("hMatchDzDist",";dZ [cm]",500,-100,100);
514 fOutput->Add(fHMatchDz);
515 fHMatchEp = new TH1F("hMatchEpDist",";E/p",100,0,10);
516 fOutput->Add(fHMatchEp);
518 // histograms for pion candidates
520 fHPionEtaPhi = new TH2F("hPionEtaPhi","",100,-0.8,0.8,100*nsm,phimin,phimax);
521 fHPionEtaPhi->SetXTitle("#eta_{#gamma#gamma}");
522 fHPionEtaPhi->SetYTitle("#varphi_{#gamma#gamma}");
523 fOutput->Add(fHPionEtaPhi);
524 fHPionMggPt = new TH2F("hPionMggPt","",1000,0,2,100,0,20.0);
525 fHPionMggPt->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
526 fHPionMggPt->SetYTitle("p_{T}^{#gamma#gamma} [GeV/c]");
527 fOutput->Add(fHPionMggPt);
528 fHPionMggAsym = new TH2F("hPionMggAsym","",1000,0,2,100,0,1);
529 fHPionMggAsym->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
530 fHPionMggAsym->SetYTitle("Z_{#gamma#gamma} [GeV]");
531 fOutput->Add(fHPionMggAsym);
532 fHPionMggDgg = new TH2F("hPionMggDgg","",1000,0,2,100,0,10);
533 fHPionMggDgg->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
534 fHPionMggDgg->SetYTitle("opening angle [grad]");
535 fOutput->Add(fHPionMggDgg);
536 const Int_t nbins = 20;
537 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};
538 fPtRanges = new TAxis(nbins-1,xbins);
539 for (Int_t i = 0; i<=nbins; ++i) {
540 fHPionInvMasses[i] = new TH1F(Form("hPionInvMass%d",i),"",1000,0,2);
541 fHPionInvMasses[i]->SetXTitle("M_{#gamma#gamma} [GeV/c^{2}]");
543 fHPionInvMasses[i]->SetTitle(Form("0 < p_{T}^{#gamma#gamma} <%.1f",xbins[0]));
545 fHPionInvMasses[i]->SetTitle(Form("p_{T}^{#gamma#gamma} > 50"));
547 fHPionInvMasses[i]->SetTitle(Form("%.1f < p_{T}^{#gamma#gamma} <%.1f",xbins[i-1],xbins[i]));
548 fOutput->Add(fHPionInvMasses[i]);
552 PostData(1, fOutput);
555 //________________________________________________________________________
556 void AliAnalysisTaskEMCALPi0PbPb::UserExec(Option_t *)
558 // Called for each event.
563 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
564 fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
565 UInt_t offtrigger = 0;
567 am->LoadBranch("AliESDRun.");
568 am->LoadBranch("AliESDHeader.");
569 UInt_t mask1 = fEsdEv->GetESDRun()->GetDetectorsInDAQ();
570 UInt_t mask2 = fEsdEv->GetESDRun()->GetDetectorsInReco();
571 Bool_t desc1 = (mask1 >> 18) & 0x1;
572 Bool_t desc2 = (mask2 >> 18) & 0x1;
573 if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
574 AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
575 mask1, fEsdEv->GetESDRun()->GetDetectorsInReco(),
576 mask2, fEsdEv->GetESDRun()->GetDetectorsInDAQ()));
579 offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
581 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
583 AliFatal("Neither ESD nor AOD event found");
586 am->LoadBranch("header");
587 offtrigger = fAodEv->GetHeader()->GetOfflineTrigger();
589 if (!fMcMode && (offtrigger & AliVEvent::kFastOnly)) {
590 AliWarning(Form("EMCAL not in fast only partition"));
593 if (fDoTrMatGeom && !AliGeomManager::GetGeometry()) { // get geometry
594 AliWarning("Accessing geometry from OCDB, this is not very efficient!");
595 AliCDBManager *cdb = AliCDBManager::Instance();
596 if (!cdb->IsDefaultStorageSet())
597 cdb->SetDefaultStorage("raw://");
598 Int_t runno = InputEvent()->GetRunNumber();
599 if (runno != cdb->GetRun())
601 AliGeomManager::LoadGeometry();
604 if (!AliGeomManager::GetGeometry()&&!fIsGeoMatsSet) { // set misalignment matrices (stored in first event)
605 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
606 for (Int_t i=0; i<nsm; ++i) {
607 const TGeoHMatrix *geom = 0;
609 geom = fEsdEv->GetESDRun()->GetEMCALMatrix(i);
611 geom = fAodEv->GetHeader()->GetEMCALMatrix(i);
615 fGeom->SetMisalMatrix(geom,i);
617 fIsGeoMatsSet = kTRUE;
620 if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
622 const AliESDRun *erun = fEsdEv->GetESDRun();
623 AliMagF *field = AliMagF::CreateFieldMap(erun->GetCurrentL3(),
624 erun->GetCurrentDip(),
627 erun->GetBeamEnergy(),
628 erun->GetBeamType());
629 TGeoGlobalMagField::Instance()->SetField(field);
631 Double_t pol = -1; //polarity
632 Double_t be = -1; //beam energy
633 AliMagF::BeamType_t btype = AliMagF::kBeamTypepp;
634 Int_t runno = fAodEv->GetRunNumber();
635 if (runno>=136851 && runno<138275) {
638 btype = AliMagF::kBeamTypeAA;
639 } else if (runno>=138275 && runno<=139517) {
642 btype = AliMagF::kBeamTypeAA;
644 AliError(Form("Do not know the bfield parameters for run %d! Using defaults!!!", runno));
646 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", pol, pol, AliMagF::k5kG, btype, be));
654 AliESDHeader *h = dynamic_cast<AliESDHeader*>(InputEvent()->GetHeader());
656 trgclasses = fEsdEv->GetFiredTriggerClasses();
658 AliAODHeader *h2 = dynamic_cast<AliAODHeader*>(InputEvent()->GetHeader());
660 trgclasses = h2->GetFiredTriggerClasses();
662 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
663 const char *name = fTrClassNamesArr->At(i)->GetName();
664 if (trgclasses.Contains(name))
665 fHTclsBeforeCuts->Fill(1+i);
668 if (fDoPSel && offtrigger==0)
673 const AliCentrality *centP = InputEvent()->GetCentrality();
674 Double_t cent = centP->GetCentralityPercentileUnchecked(fCentVar);
676 if (cent<fCentFrom||cent>fCentTo)
682 if (centP->GetQuality()>0)
686 fHCentQual->Fill(cent);
690 am->LoadBranch("PrimaryVertex.");
691 am->LoadBranch("SPDVertex.");
692 am->LoadBranch("TPCVertex.");
694 fAodEv = dynamic_cast<AliAODEvent*>(InputEvent());
695 am->LoadBranch("vertices");
699 const AliVVertex *vertex = InputEvent()->GetPrimaryVertex();
703 fHVertexZ->Fill(vertex->GetZ());
705 if(vertex->GetZ()<fVtxZMin||vertex->GetZ()>fVtxZMax)
709 fHVertexZ2->Fill(vertex->GetZ());
711 // count number of accepted events
714 for (Int_t i = 0; i<fTrClassNamesArr->GetEntries(); ++i) {
715 const char *name = fTrClassNamesArr->At(i)->GetName();
716 if (trgclasses.Contains(name))
717 fHTclsAfterCuts->Fill(1+i);
720 fRecPoints = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
721 fDigits = 0; // will be set if fClusName is given and AliAnalysisTaskEMCALClusterizeFast is used
722 fEsdClusters = 0; // will be set if ESD input used and if fRecPoints are not set or if clusters are attached
723 fEsdCells = 0; // will be set if ESD input used
724 fAodClusters = 0; // will be set if AOD input used and if fRecPoints are not set or if clusters are attached
725 // or if fClusName is given and AliAnalysisTaskEMCALClusterizeFast in AOD output mode
726 fAodCells = 0; // will be set if AOD input used
728 // deal with special output from AliAnalysisTaskEMCALClusterizeFast first
729 Bool_t clusattached = 0;
730 Bool_t recalibrated = 0;
731 if (1 && !fClusName.IsNull()) {
732 AliAnalysisTaskEMCALClusterizeFast *cltask = 0;
733 TObjArray *ts = am->GetTasks();
734 cltask = dynamic_cast<AliAnalysisTaskEMCALClusterizeFast*>(ts->FindObject(fClusName));
735 if (cltask && cltask->GetClusters()) {
736 fRecPoints = cltask->GetClusters();
737 fDigits = cltask->GetDigits();
738 clusattached = cltask->GetAttachClusters();
739 if (cltask->GetCalibData()!=0)
740 recalibrated = kTRUE;
743 if (1 && !fClusName.IsNull()) {
746 l = AODEvent()->GetList();
748 l = fAodEv->GetList();
750 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject(fClusName));
754 if (fEsdEv) { // ESD input mode
755 if (1 && (!fRecPoints||clusattached)) {
757 am->LoadBranch("CaloClusters");
758 TList *l = fEsdEv->GetList();
760 fEsdClusters = dynamic_cast<TClonesArray*>(l->FindObject("CaloClusters"));
765 am->LoadBranch("EMCALCells.");
766 fEsdCells = fEsdEv->GetEMCALCells();
768 } else if (fAodEv) { // AOD input mode
769 if (1 && (!fAodClusters || clusattached)) {
771 am->LoadBranch("caloClusters");
772 TList *l = fAodEv->GetList();
774 fAodClusters = dynamic_cast<TClonesArray*>(l->FindObject("caloClusters"));
779 am->LoadBranch("emcalCells");
780 fAodCells = fAodEv->GetEMCALCells();
783 AliFatal("Impossible to not have either pointer to ESD or AOD event");
787 AliDebug(2,Form("fRecPoints set: %p", fRecPoints));
788 AliDebug(2,Form("fDigits set: %p", fDigits));
789 AliDebug(2,Form("fEsdClusters set: %p", fEsdClusters));
790 AliDebug(2,Form("fEsdCells set: %p", fEsdCells));
791 AliDebug(2,Form("fAodClusters set: %p", fAodClusters));
792 AliDebug(2,Form("fAodCells set: %p", fAodCells));
796 ClusterAfterburner();
815 fSelPrimTracks->Clear();
824 PostData(1, fOutput);
827 //________________________________________________________________________
828 void AliAnalysisTaskEMCALPi0PbPb::Terminate(Option_t *)
830 // Terminate called at the end of analysis.
833 TFile *f = OpenFile(1);
834 TDirectory::TContext context(f);
839 AliInfo(Form("%s: Accepted %lld events ", GetName(), fNEvs));
842 //________________________________________________________________________
843 void AliAnalysisTaskEMCALPi0PbPb::CalcCaloTriggers()
845 // Calculate triggers
848 return; // information not available in AOD
855 AliVCaloCells *cells = fEsdCells;
861 Int_t ncells = cells->GetNumberOfCells();
865 if (ncells>fNAmpInTrigger) {
866 delete [] fAmpInTrigger;
867 fAmpInTrigger = new Float_t[ncells];
868 fNAmpInTrigger = ncells;
870 for (Int_t i=0;i<ncells;++i)
871 fAmpInTrigger[i] = 0;
873 std::map<Short_t,Short_t> map;
874 for (Short_t pos=0;pos<ncells;++pos) {
875 Short_t id = cells->GetCellNumber(pos);
879 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
880 am->LoadBranch("EMCALTrigger.");
882 AliESDCaloTrigger *triggers = fEsdEv->GetCaloTrigger("EMCAL");
885 if (triggers->GetEntries()<=0)
890 while (triggers->Next()) {
891 Int_t gCol=0, gRow=0, ntimes=0;
892 triggers->GetPosition(gCol,gRow);
893 triggers->GetNL0Times(ntimes);
897 triggers->GetAmplitude(amp);
899 fGeom->GetAbsFastORIndexFromPositionInEMCAL(gCol,gRow,find);
902 Int_t cidx[4] = {-1};
903 Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
907 triggers->GetL0Times(trgtimes);
908 Int_t mintime = trgtimes[0];
909 Int_t maxtime = trgtimes[0];
910 Bool_t trigInTimeWindow = 0;
911 for (Int_t i=0;i<ntimes;++i) {
912 if (trgtimes[i]<mintime)
913 mintime = trgtimes[i];
914 if (maxtime<trgtimes[i])
915 maxtime = trgtimes[i];
916 if ((fMinL0Time<=trgtimes[i]) && (fMaxL0Time>=trgtimes[i]))
917 trigInTimeWindow = 1;
920 Double_t tenergy = 0;
923 for (Int_t i=0;i<3;++i) {
925 std::map<Short_t,Short_t>::iterator it = map.find(cidx[i]);
930 if (trigInTimeWindow)
931 fAmpInTrigger[pos] = amp;
932 Float_t eta=-1, phi=-1;
933 fGeom->EtaPhiFromIndex(cidx[i],eta,phi);
934 Double_t en= cells->GetAmplitude(pos);
946 AliStaTrigger *trignew = static_cast<AliStaTrigger*>(fTriggers->New(ntrigs++));
947 trignew->fE = tenergy;
948 trignew->fEta = teta;
949 trignew->fPhi = tphi;
951 trignew->fMinTime = mintime;
952 trignew->fMaxTime = maxtime;
956 //________________________________________________________________________
957 void AliAnalysisTaskEMCALPi0PbPb::CalcClusterProps()
959 // Calculate cluster properties
966 AliVCaloCells *cells = fEsdCells;
972 TObjArray *clusters = fEsdClusters;
974 clusters = fAodClusters;
978 Int_t ncells = cells->GetNumberOfCells();
979 Int_t nclus = clusters->GetEntries();
980 Int_t ntrks = fSelTracks->GetEntries();
981 Bool_t btracks[6][ntrks];
982 memset(btracks,0,sizeof(btracks));
984 std::map<Short_t,Short_t> map;
985 for (Short_t pos=0;pos<ncells;++pos) {
986 Short_t id = cells->GetCellNumber(pos);
990 TObjArray filtMcParts;
992 Int_t nmc = fMcParts->GetEntries();
993 for (Int_t i=0; i<nmc; ++i) {
994 AliStaPart *pa = static_cast<AliStaPart*>(fMcParts->At(i));
1000 for(Int_t i=0, ncl=0; i<nclus; ++i) {
1001 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1005 if (!clus->IsEMCAL())
1007 if (clus->E()<fMinE)
1010 Float_t clsPos[3] = {0};
1011 clus->GetPosition(clsPos);
1012 TVector3 clsVec(clsPos);
1013 Double_t vertex[3] = {0};
1014 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1015 TLorentzVector clusterVec;
1016 clus->GetMomentum(clusterVec,vertex);
1017 Double_t clsEta = clusterVec.Eta();
1019 AliStaCluster *cl = static_cast<AliStaCluster*>(fClusters->New(ncl++));
1021 cl->fR = clsVec.Perp();
1022 cl->fEta = clsVec.Eta();
1023 cl->fPhi = clsVec.Phi();
1024 cl->fN = clus->GetNCells();
1025 cl->fN1 = GetNCells(clus,0.100);
1026 cl->fN3 = GetNCells(clus,0.300);
1028 Double_t emax = GetMaxCellEnergy(clus, id);
1031 cl->fTmax = cells->GetCellTime(id);
1032 if (clus->GetDistanceToBadChannel()<10000)
1033 cl->fDbc = clus->GetDistanceToBadChannel();
1034 if (!TMath::IsNaN(clus->GetDispersion()))
1035 cl->fDisp = clus->GetDispersion();
1036 if (!TMath::IsNaN(clus->GetM20()))
1037 cl->fM20 = clus->GetM20();
1038 if (!TMath::IsNaN(clus->GetM02()))
1039 cl->fM02 = clus->GetM02();
1040 Double_t maxAxis = 0, minAxis = 0;
1041 GetSigma(clus,maxAxis,minAxis);
1042 clus->SetTOF(maxAxis); // store sigma in TOF
1044 Double_t clusterEcc = 0;
1046 clusterEcc = TMath::Sqrt(1.0 - minAxis*minAxis/(maxAxis*maxAxis));
1047 clus->SetChi2(clusterEcc); // store ecc in chi2
1048 cl->fEcc = clusterEcc;
1049 cl->fTrIso = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist);
1050 cl->fTrIso1 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 1);
1051 cl->fTrIso2 = GetTrackIsolation(clusterVec.Eta(),clusterVec.Phi(),fIsoDist, 2);
1052 cl->fCeCore = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),0.05);
1053 cl->fCeIso = GetCellIsolation(clsVec.Eta(),clsVec.Phi(),fIsoDist);
1055 if (fAmpInTrigger) { // fill trigger info if present
1056 Double_t trigpen = 0;
1057 Double_t trignen = 0;
1058 for(Int_t j=0; j<cl->fN; ++j) {
1059 Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
1061 std::map<Short_t,Short_t>::iterator it = map.find(cid);
1066 if (fAmpInTrigger[pos]>0)
1067 trigpen += cells->GetAmplitude(pos);
1068 else if (fAmpInTrigger[pos]<0)
1069 trignen += cells->GetAmplitude(pos);
1073 cl->fTrigE = trigpen;
1077 cl->fTrigMaskE = trignen;
1080 cl->fIsShared = IsShared(clus);
1083 Double_t mind2 = 1e10;
1084 for(Int_t j = 0; j<ntrks; ++j) {
1085 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
1089 if (TMath::Abs(clsEta-track->Eta())>0.5)
1092 TVector3 vec(clsPos);
1093 Int_t index = (Int_t)(vec.Phi()*TMath::RadToDeg()/20);
1094 if (btracks[index-4][j]) {
1098 Float_t tmpR=-1, tmpZ=-1;
1099 if (!fDoTrMatGeom) {
1100 AliExternalTrackParam *tParam = 0;
1102 AliESDtrack *esdTrack = static_cast<AliESDtrack*>(track);
1103 tParam = new AliExternalTrackParam(*esdTrack->GetTPCInnerParam());
1105 tParam = new AliExternalTrackParam(track);
1107 Double_t bfield[3] = {0};
1108 track->GetBxByBz(bfield);
1109 Double_t alpha = (index+0.5)*20*TMath::DegToRad();
1110 vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
1111 tParam->Rotate(alpha); //Rotate the track to the same local extrapolation system
1112 Bool_t ret = tParam->PropagateToBxByBz(vec.X(), bfield);
1114 btracks[index-4][j]=1;
1118 Double_t trkPos[3] = {0};
1119 tParam->GetXYZ(trkPos); //Get the extrapolated global position
1120 tmpR = TMath::Sqrt( TMath::Power(clsPos[0]-trkPos[0],2) +
1121 TMath::Power(clsPos[1]-trkPos[1],2) +
1122 TMath::Power(clsPos[2]-trkPos[2],2) );
1123 tmpZ = clsPos[2]-trkPos[2];
1126 if (TMath::Abs(clsEta-track->Eta())>fIsoDist)
1128 AliExternalTrackParam tParam(track);
1129 if (!fReco->ExtrapolateTrackToCluster(&tParam, clus, tmpR, tmpZ))
1137 cl->fTrDr = TMath::Sqrt(tmpR*tmpR-tmpZ*tmpZ);
1138 cl->fTrEp = clus->E()/track->P();
1143 if (cl->fIsTrackM) {
1144 fHMatchDr->Fill(cl->fTrDr);
1145 fHMatchDz->Fill(cl->fTrDz);
1146 fHMatchEp->Fill(cl->fTrEp);
1151 Int_t nmc = filtMcParts.GetEntries();
1152 Double_t diffR2 = 1e9;
1153 AliStaPart *msta = 0;
1154 for (Int_t j=0; j<nmc; ++j) {
1155 AliStaPart *pa = static_cast<AliStaPart*>(filtMcParts.At(j));
1156 Double_t t1=clsVec.Eta()-pa->fVEta;
1157 Double_t t2=TVector2::Phi_mpi_pi(clsVec.Phi()-pa->fVPhi);
1158 Double_t tmp = t1*t1+t2*t2;
1164 if (diffR2<10 && msta!=0) {
1165 cl->fMcLabel = msta->fLab;
1170 if (fDigits && fEmbedMode) {
1171 for(Int_t j=0; j<cl->fN; ++j) {
1172 Short_t cid = TMath::Abs(clus->GetCellAbsId(j));
1174 std::map<Short_t,Short_t>::iterator it = map.find(cid);
1179 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigits->At(pos));
1182 if (digit->GetId() != cid) {
1183 AliError(Form("Ids should be equal: %d %d", cid, digit->GetId()));
1186 if (digit->GetType()<-1) {
1187 cl->fEmbE += digit->GetChi2();
1194 //________________________________________________________________________
1195 void AliAnalysisTaskEMCALPi0PbPb::CalcPrimTracks()
1197 // Calculate track properties.
1199 fSelPrimTracks->Clear();
1201 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1202 TClonesArray *tracks = 0;
1204 am->LoadBranch("Tracks");
1205 TList *l = fEsdEv->GetList();
1206 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1208 am->LoadBranch("tracks");
1209 TList *l = fAodEv->GetList();
1210 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1216 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1217 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad()-fIsoDist*1.25;
1218 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad()+fIsoDist*1.25;
1221 fSelPrimTracks->SetOwner(kTRUE);
1222 am->LoadBranch("PrimaryVertex.");
1223 const AliESDVertex *vtx = fEsdEv->GetPrimaryVertexSPD();
1224 am->LoadBranch("SPDVertex.");
1225 const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
1226 am->LoadBranch("Tracks");
1227 const Int_t Ntracks = tracks->GetEntries();
1228 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1229 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1231 AliWarning(Form("Could not receive track %d\n", iTracks));
1234 if (fTrCuts && !fTrCuts->IsSelected(track))
1236 Double_t eta = track->Eta();
1239 if (track->Phi()<phimin||track->Phi()>phimax)
1242 AliESDtrack copyt(*track);
1244 copyt.GetBxByBz(bfield);
1245 AliExternalTrackParam tParam;
1246 Bool_t relate = copyt.RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&tParam);
1249 copyt.Set(tParam.GetX(),tParam.GetAlpha(),tParam.GetParameter(),tParam.GetCovariance());
1251 Double_t p[3] = { 0. };
1253 Double_t pos[3] = { 0. };
1255 Double_t covTr[21] = { 0. };
1256 copyt.GetCovarianceXYZPxPyPz(covTr);
1257 Double_t pid[10] = { 0. };
1258 copyt.GetESDpid(pid);
1259 AliAODTrack *aTrack = new AliAODTrack(copyt.GetID(),
1266 (Short_t)copyt.GetSign(),
1267 copyt.GetITSClusterMap(),
1269 0,/*fPrimaryVertex,*/
1270 kTRUE, // check if this is right
1271 vtx->UsesTrack(copyt.GetID()));
1272 aTrack->SetTPCClusterMap(copyt.GetTPCClusterMap());
1273 aTrack->SetTPCSharedMap (copyt.GetTPCSharedMap());
1274 Float_t ndf = copyt.GetTPCNcls() + 1 - 5 ;
1276 aTrack->SetChi2perNDF(copyt.GetTPCchi2()/ndf);
1278 aTrack->SetChi2perNDF(-1);
1279 aTrack->SetFlags(copyt.GetStatus());
1280 aTrack->SetTPCPointsF(copyt.GetTPCNclsF());
1281 fSelPrimTracks->Add(aTrack);
1284 Int_t ntracks = tracks->GetEntries();
1285 for (Int_t i=0; i<ntracks; ++i) {
1286 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1289 Double_t eta = track->Eta();
1292 if (track->Phi()<phimin||track->Phi()>phimax)
1294 if(track->GetTPCNcls()<fMinNClusPerTr)
1296 //todo: Learn how to set/filter AODs for prim/sec tracks
1297 fSelPrimTracks->Add(track);
1302 //________________________________________________________________________
1303 void AliAnalysisTaskEMCALPi0PbPb::CalcTracks()
1305 // Calculate track properties (including secondaries).
1307 fSelTracks->Clear();
1309 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1310 TClonesArray *tracks = 0;
1312 am->LoadBranch("Tracks");
1313 TList *l = fEsdEv->GetList();
1314 tracks = dynamic_cast<TClonesArray*>(l->FindObject("Tracks"));
1316 am->LoadBranch("tracks");
1317 TList *l = fAodEv->GetList();
1318 tracks = dynamic_cast<TClonesArray*>(l->FindObject("tracks"));
1325 const Int_t Ntracks = tracks->GetEntries();
1326 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1327 AliESDtrack *track = static_cast<AliESDtrack*>(tracks->At(iTracks));
1329 AliWarning(Form("Could not receive track %d\n", iTracks));
1332 if (fTrCuts && !fTrCuts->IsSelected(track))
1334 Double_t eta = track->Eta();
1337 fSelTracks->Add(track);
1340 Int_t ntracks = tracks->GetEntries();
1341 for (Int_t i=0; i<ntracks; ++i) {
1342 AliAODTrack *track = static_cast<AliAODTrack*>(tracks->At(i));
1345 Double_t eta = track->Eta();
1348 if(track->GetTPCNcls()<fMinNClusPerTr)
1351 if (0 && (track->Pt()>=0.6) && (track->PxAtDCA()==-999)) { // compute position on EMCAL
1352 AliExternalTrackParam tParam(track);
1353 if (AliTrackerBase::PropagateTrackToBxByBz(&tParam, 438, 0.139, 1, kTRUE)) {
1355 tParam.GetXYZ(trkPos);
1356 track->SetPxPyPzAtDCA(trkPos[0],trkPos[1],trkPos[2]);
1359 fSelTracks->Add(track);
1364 //________________________________________________________________________
1365 void AliAnalysisTaskEMCALPi0PbPb::ClusterAfterburner()
1367 // Run custer reconstruction afterburner.
1369 AliVCaloCells *cells = fEsdCells;
1376 Int_t ncells = cells->GetNumberOfCells();
1380 Double_t cellMeanE = 0, cellSigE = 0;
1381 for (Int_t i = 0; i<ncells; ++i) {
1382 Double_t cellE = cells->GetAmplitude(i);
1384 cellSigE += cellE*cellE;
1386 cellMeanE /= ncells;
1388 cellSigE -= cellMeanE*cellMeanE;
1391 cellSigE = TMath::Sqrt(cellSigE / ncells);
1393 Double_t subE = cellMeanE - 7*cellSigE;
1397 for (Int_t i = 0; i<ncells; ++i) {
1399 Double_t amp=0,time=0;
1400 if (!cells->GetCell(i, id, amp, time))
1405 cells->SetCell(i, id, amp, time);
1408 TObjArray *clusters = fEsdClusters;
1410 clusters = fAodClusters;
1414 Int_t nclus = clusters->GetEntries();
1415 for (Int_t i = 0; i<nclus; ++i) {
1416 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1417 if (!clus->IsEMCAL())
1419 Int_t nc = clus->GetNCells();
1421 UShort_t ids[100] = {0};
1422 Double_t fra[100] = {0};
1423 for (Int_t j = 0; j<nc; ++j) {
1424 Short_t id = TMath::Abs(clus->GetCellAbsId(j));
1425 Double_t cen = cells->GetCellAmplitude(id);
1433 clusters->RemoveAt(i);
1437 for (Int_t j = 0; j<nc; ++j) {
1438 Short_t id = ids[j];
1439 Double_t cen = cells->GetCellAmplitude(id);
1443 AliAODCaloCluster *aodclus = dynamic_cast<AliAODCaloCluster*>(clus);
1446 aodclus->SetNCells(nc);
1447 aodclus->SetCellsAmplitudeFraction(fra);
1448 aodclus->SetCellsAbsId(ids);
1451 clusters->Compress();
1454 //________________________________________________________________________
1455 void AliAnalysisTaskEMCALPi0PbPb::FillCellHists()
1457 // Fill histograms related to cell properties.
1459 AliVCaloCells *cells = fEsdCells;
1466 Int_t cellModCount[12] = {0};
1467 Double_t cellMaxE = 0;
1468 Double_t cellMeanE = 0;
1469 Int_t ncells = cells->GetNumberOfCells();
1473 Int_t nsm = fGeom->GetEMCGeometry()->GetNumberOfSuperModules();
1475 for (Int_t i = 0; i<ncells; ++i) {
1476 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1477 Double_t cellE = cells->GetAmplitude(i);
1478 fHCellE->Fill(cellE);
1483 Int_t iSM=-1, iTower=-1, nIphi=-1, nIeta=-1;
1484 Bool_t ret = fGeom->GetCellIndex(absID, iSM, iTower, nIphi, nIeta);
1486 AliError(Form("Could not get cell index for %d", absID));
1489 ++cellModCount[iSM];
1490 Int_t iPhi=-1, iEta=-1;
1491 fGeom->GetCellPhiEtaIndexInSModule(iSM, iTower, nIphi, nIeta, iPhi, iEta);
1492 fHColuRow[iSM]->Fill(iEta,iPhi,1);
1493 fHColuRowE[iSM]->Fill(iEta,iPhi,cellE);
1494 fHCellFreqNoCut[iSM]->Fill(absID);
1495 if (cellE > 0.1) fHCellFreqCut100M[iSM]->Fill(absID);
1496 if (cellE > 0.3) fHCellFreqCut300M[iSM]->Fill(absID);
1497 if (fHCellCheckE && fHCellCheckE[absID])
1498 fHCellCheckE[absID]->Fill(cellE);
1499 fHCellFreqE[iSM]->Fill(absID, cellE);
1501 fHCellH->Fill(cellMaxE);
1502 cellMeanE /= ncells;
1503 fHCellM->Fill(cellMeanE);
1504 fHCellM2->Fill(cellMeanE*ncells/24/48/nsm);
1505 for (Int_t i=0; i<nsm; ++i)
1506 fHCellMult[i]->Fill(cellModCount[i]);
1509 //________________________________________________________________________
1510 void AliAnalysisTaskEMCALPi0PbPb::FillClusHists()
1512 // Fill histograms related to cluster properties.
1514 TObjArray *clusters = fEsdClusters;
1516 clusters = fAodClusters;
1520 Double_t vertex[3] = {0};
1521 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1523 Int_t nclus = clusters->GetEntries();
1524 for(Int_t i = 0; i<nclus; ++i) {
1525 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(i));
1528 if (!clus->IsEMCAL())
1530 TLorentzVector clusterVec;
1531 clus->GetMomentum(clusterVec,vertex);
1532 Double_t maxAxis = clus->GetTOF(); //sigma
1533 Double_t clusterEcc = clus->Chi2(); //eccentricity
1534 fHClustEccentricity->Fill(clusterEcc);
1535 fHClustEtaPhi->Fill(clusterVec.Eta(),clusterVec.Phi());
1536 fHClustEnergyPt->Fill(clusterVec.E(),clusterVec.Pt());
1537 fHClustEnergySigma->Fill(clus->E()*maxAxis,clus->E());
1538 fHClustSigmaSigma->Fill(max(clus->GetM02(),clus->GetM20()),clus->E()*maxAxis);
1539 fHClustNCellEnergyRatio->Fill(clus->GetNCells(),GetMaxCellEnergy(clus)/clus->E());
1541 fHClustEnergyNCell->Fill(clus->E(),clus->GetNCells());
1545 //________________________________________________________________________
1546 void AliAnalysisTaskEMCALPi0PbPb::CalcMcInfo()
1548 // Get Mc truth particle information.
1555 AliEMCALEMCGeometry *emc = fGeom->GetEMCGeometry();
1556 Double_t etamin = -0.7;
1557 Double_t etamax = +0.7;
1558 Double_t phimin = emc->GetArm1PhiMin()*TMath::DegToRad();
1559 Double_t phimax = emc->GetArm1PhiMax()*TMath::DegToRad();
1562 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1563 am->LoadBranch(AliAODMCParticle::StdBranchName());
1564 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAodEv->FindListObject(AliAODMCParticle::StdBranchName()));
1568 Int_t nents = tca->GetEntries();
1569 for(int it=0; it<nents; ++it) {
1570 AliAODMCParticle *part = static_cast<AliAODMCParticle*>(tca->At(it));
1573 // pion or eta meson or direct photon
1574 if(part->GetPdgCode() == 111) {
1575 } else if(part->GetPdgCode() == 221) {
1576 } else if(part->GetPdgCode() == 22 ) {
1581 Double_t dR = TMath::Sqrt((part->Xv()*part->Xv())+(part->Yv()*part->Yv()));
1586 Double_t pt = part->Pt() ;
1589 Double_t eta = part->Eta();
1590 if (eta<etamin||eta>etamax)
1592 Double_t phi = part->Phi();
1593 if (phi<phimin||phi>phimax)
1596 ProcessDaughters(part, it, tca);
1601 AliMCEvent *mcEvent = MCEvent();
1605 const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
1609 mcEvent->PreReadAll();
1611 Int_t nTracks = mcEvent->GetNumberOfPrimaries();
1612 for (Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
1613 AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1617 // pion or eta meson or direct photon
1618 if(mcP->PdgCode() == 111) {
1619 } else if(mcP->PdgCode() == 221) {
1620 } else if(mcP->PdgCode() == 22 ) {
1625 if(mcP->GetMother()>=0 && mcP->GetMother()<nTracks)
1627 Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1628 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1633 Double_t pt = mcP->Pt() ;
1636 Double_t eta = mcP->Eta();
1637 if (eta<etamin||eta>etamax)
1639 Double_t phi = mcP->Phi();
1640 if (phi<phimin||phi>phimax)
1643 ProcessDaughters(mcP, iTrack, mcEvent);
1647 //________________________________________________________________________
1648 void AliAnalysisTaskEMCALPi0PbPb::FillNtuple()
1655 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
1657 fHeader->fRun = fAodEv->GetRunNumber();
1658 fHeader->fOrbit = fAodEv->GetHeader()->GetOrbitNumber();
1659 fHeader->fPeriod = fAodEv->GetHeader()->GetPeriodNumber();
1660 fHeader->fBx = fAodEv->GetHeader()->GetBunchCrossNumber();
1661 fHeader->fL0 = fAodEv->GetHeader()->GetL0TriggerInputs();
1662 fHeader->fL1 = fAodEv->GetHeader()->GetL1TriggerInputs();
1663 fHeader->fL2 = fAodEv->GetHeader()->GetL2TriggerInputs();
1664 fHeader->fTrClassMask = fAodEv->GetHeader()->GetTriggerMask();
1665 fHeader->fTrCluster = fAodEv->GetHeader()->GetTriggerCluster();
1666 fHeader->fOffTriggers = fAodEv->GetHeader()->GetOfflineTrigger();
1667 fHeader->fFiredTriggers = fAodEv->GetHeader()->GetFiredTriggerClasses();
1669 fHeader->fRun = fEsdEv->GetRunNumber();
1670 fHeader->fOrbit = fEsdEv->GetHeader()->GetOrbitNumber();
1671 fHeader->fPeriod = fEsdEv->GetESDRun()->GetPeriodNumber();
1672 fHeader->fBx = fEsdEv->GetHeader()->GetBunchCrossNumber();
1673 fHeader->fL0 = fEsdEv->GetHeader()->GetL0TriggerInputs();
1674 fHeader->fL1 = fEsdEv->GetHeader()->GetL1TriggerInputs();
1675 fHeader->fL2 = fEsdEv->GetHeader()->GetL2TriggerInputs();
1676 fHeader->fTrClassMask = fEsdEv->GetHeader()->GetTriggerMask();
1677 fHeader->fTrCluster = fEsdEv->GetHeader()->GetTriggerCluster();
1678 fHeader->fOffTriggers = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
1679 fHeader->fFiredTriggers = fEsdEv->GetFiredTriggerClasses();
1680 Float_t v0CorrR = 0;
1681 fHeader->fV0 = AliESDUtils::GetCorrV0(fEsdEv,v0CorrR);
1682 const AliMultiplicity *mult = fEsdEv->GetMultiplicity();
1684 fHeader->fCl1 = mult->GetNumberOfITSClusters(1);
1685 fHeader->fTr = AliESDtrackCuts::GetReferenceMultiplicity(fEsdEv,1);
1687 AliCentrality *cent = InputEvent()->GetCentrality();
1688 fHeader->fV0Cent = cent->GetCentralityPercentileUnchecked("V0M");
1689 fHeader->fCl1Cent = cent->GetCentralityPercentileUnchecked("CL1");
1690 fHeader->fTrCent = cent->GetCentralityPercentileUnchecked("TRK");
1691 fHeader->fCqual = cent->GetQuality();
1693 AliEventplane *ep = InputEvent()->GetEventplane();
1695 if (ep->GetQVector())
1696 fHeader->fPsi = ep->GetQVector()->Phi()/2. ;
1699 if (ep->GetQsub1()&&ep->GetQsub2())
1700 fHeader->fPsiRes = ep->GetQsub1()->Phi()/2.-ep->GetQsub2()->Phi()/2.;
1702 fHeader->fPsiRes = 0;
1706 TString trgclasses(fHeader->fFiredTriggers);
1707 for (Int_t j = 0; j<fTrClassNamesArr->GetEntries(); ++j) {
1708 const char *name = fTrClassNamesArr->At(j)->GetName();
1709 if (trgclasses.Contains(name))
1710 val += TMath::Power(2,j);
1712 fHeader->fTcls = (UInt_t)val;
1714 fHeader->fNSelTr = fSelTracks->GetEntries();
1715 fHeader->fNSelPrimTr = fSelPrimTracks->GetEntries();
1716 fHeader->fNSelPrimTr1 = 0;
1717 fHeader->fNSelPrimTr2 = 0;
1718 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++){
1719 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1721 ++fHeader->fNSelPrimTr1;
1723 ++fHeader->fNSelPrimTr2;
1726 fHeader->fNCells = 0;
1727 fHeader->fNCells1 = 0;
1728 fHeader->fNCells2 = 0;
1729 fHeader->fNCells5 = 0;
1730 fHeader->fMaxCellE = 0;
1732 AliVCaloCells *cells = fEsdCells;
1737 Int_t ncells = cells->GetNumberOfCells();
1738 for(Int_t j=0; j<ncells; ++j) {
1739 Double_t cellen = cells->GetAmplitude(j);
1741 ++fHeader->fNCells1;
1743 ++fHeader->fNCells2;
1745 ++fHeader->fNCells5;
1746 if (cellen>fHeader->fMaxCellE)
1747 fHeader->fMaxCellE = cellen;
1749 fHeader->fNCells = ncells;
1752 fHeader->fNClus = 0;
1753 fHeader->fNClus1 = 0;
1754 fHeader->fNClus2 = 0;
1755 fHeader->fNClus5 = 0;
1756 fHeader->fMaxClusE = 0;
1758 TObjArray *clusters = fEsdClusters;
1760 clusters = fAodClusters;
1763 Int_t nclus = clusters->GetEntries();
1764 for(Int_t j=0; j<nclus; ++j) {
1765 AliVCluster *clus = static_cast<AliVCluster*>(clusters->At(j));
1766 if (!clus->IsEMCAL())
1768 Double_t clusen = clus->E();
1775 if (clusen>fHeader->fMaxClusE)
1776 fHeader->fMaxClusE = clusen;
1778 fHeader->fNClus = nclus;
1782 am->LoadBranch("vertices");
1783 AliAODVertex *pv = fAodEv->GetPrimaryVertex();
1784 FillVertex(fPrimVert, pv);
1785 AliAODVertex *sv = fAodEv->GetPrimaryVertexSPD();
1786 FillVertex(fSpdVert, sv);
1788 am->LoadBranch("PrimaryVertex.");
1789 const AliESDVertex *pv = fEsdEv->GetPrimaryVertexTracks();
1790 FillVertex(fPrimVert, pv);
1791 am->LoadBranch("SPDVertex.");
1792 const AliESDVertex *sv = fEsdEv->GetPrimaryVertexSPD();
1793 FillVertex(fSpdVert, sv);
1794 am->LoadBranch("TPCVertex.");
1795 const AliESDVertex *tv = fEsdEv->GetPrimaryVertexTPC();
1796 FillVertex(fTpcVert, tv);
1802 //________________________________________________________________________
1803 void AliAnalysisTaskEMCALPi0PbPb::FillPionHists()
1805 // Fill histograms related to pions.
1807 Double_t vertex[3] = {0};
1808 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1810 TObjArray *clusters = fEsdClusters;
1812 clusters = fAodClusters;
1815 TLorentzVector clusterVec1;
1816 TLorentzVector clusterVec2;
1817 TLorentzVector pionVec;
1819 Int_t nclus = clusters->GetEntries();
1820 for (Int_t i = 0; i<nclus; ++i) {
1821 AliVCluster *clus1 = static_cast<AliVCluster*>(clusters->At(i));
1824 if (!clus1->IsEMCAL())
1826 if (clus1->E()<fMinE)
1828 if (clus1->GetNCells()<fNminCells)
1830 if (GetMaxCellEnergy(clus1)/clus1->E()<fMinErat)
1832 if (clus1->Chi2()<fMinEcc) // eccentricity cut
1834 clus1->GetMomentum(clusterVec1,vertex);
1835 for (Int_t j = i+1; j<nclus; ++j) {
1836 AliVCluster *clus2 = static_cast<AliVCluster*>(clusters->At(j));
1839 if (!clus2->IsEMCAL())
1841 if (clus2->E()<fMinE)
1843 if (clus2->GetNCells()<fNminCells)
1845 if (GetMaxCellEnergy(clus2)/clus2->E()<fMinErat)
1847 if (clus2->Chi2()<fMinEcc) // eccentricity cut
1849 clus2->GetMomentum(clusterVec2,vertex);
1850 pionVec = clusterVec1 + clusterVec2;
1851 Double_t pionZgg = TMath::Abs(clusterVec1.E()-clusterVec2.E())/pionVec.E();
1852 Double_t pionOpeningAngle = clusterVec1.Angle(clusterVec2.Vect());
1853 if (pionZgg < fAsymMax) {
1854 fHPionEtaPhi->Fill(pionVec.Eta(),pionVec.Phi());
1855 fHPionMggPt->Fill(pionVec.M(),pionVec.Pt());
1856 fHPionMggAsym->Fill(pionVec.M(),pionZgg);
1857 fHPionMggDgg->Fill(pionVec.M(),pionOpeningAngle);
1858 Int_t bin = fPtRanges->FindBin(pionVec.Pt());
1859 fHPionInvMasses[bin]->Fill(pionVec.M());
1866 //________________________________________________________________________
1867 void AliAnalysisTaskEMCALPi0PbPb::FillMcHists()
1869 // Fill additional MC information histograms.
1874 // check if aod or esd mc mode and the fill histos
1877 //________________________________________________________________________
1878 void AliAnalysisTaskEMCALPi0PbPb::FillOtherHists()
1880 // Fill other histograms.
1882 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); ++iTracks){
1883 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1886 fHPrimTrackPt->Fill(track->Pt());
1887 fHPrimTrackEta->Fill(track->Eta());
1888 fHPrimTrackPhi->Fill(track->Phi());
1892 //________________________________________________________________________
1893 void AliAnalysisTaskEMCALPi0PbPb::FillTrackHists()
1895 // Fill track histograms.
1897 if (fSelPrimTracks) {
1898 for(int iTracks=0; iTracks < fSelPrimTracks->GetEntries(); iTracks++) {
1899 AliESDtrack *track = static_cast<AliESDtrack*>(fSelPrimTracks->At(iTracks));
1902 fHPrimTrackPt->Fill(track->Pt());
1903 fHPrimTrackEta->Fill(track->Eta());
1904 fHPrimTrackPhi->Fill(track->Phi());
1909 //__________________________________________________________________________________________________
1910 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliESDVertex *esdv)
1912 // Fill vertex from ESD vertex info.
1914 v->fVx = esdv->GetX();
1915 v->fVy = esdv->GetY();
1916 v->fVz = esdv->GetZ();
1917 v->fVc = esdv->GetNContributors();
1918 v->fDisp = esdv->GetDispersion();
1919 v->fZres = esdv->GetZRes();
1920 v->fChi2 = esdv->GetChi2();
1921 v->fSt = esdv->GetStatus();
1922 v->fIs3D = esdv->IsFromVertexer3D();
1923 v->fIsZ = esdv->IsFromVertexerZ();
1926 //__________________________________________________________________________________________________
1927 void AliAnalysisTaskEMCALPi0PbPb::FillVertex(AliStaVertex *v, const AliAODVertex *aodv)
1929 // Fill vertex from AOD vertex info.
1931 v->fVx = aodv->GetX();
1932 v->fVy = aodv->GetY();
1933 v->fVz = aodv->GetZ();
1934 v->fVc = aodv->GetNContributors();
1935 v->fChi2 = aodv->GetChi2();
1938 //________________________________________________________________________
1939 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellIsolation(Double_t cEta, Double_t cPhi, Double_t radius) const
1941 // Compute isolation based on cell content.
1943 AliVCaloCells *cells = fEsdCells;
1949 Double_t cellIsolation = 0;
1950 Double_t rad2 = radius*radius;
1951 Int_t ncells = cells->GetNumberOfCells();
1952 for (Int_t i = 0; i<ncells; ++i) {
1953 Int_t absID = TMath::Abs(cells->GetCellNumber(i));
1954 Float_t eta=-1, phi=-1;
1955 fGeom->EtaPhiFromIndex(absID,eta,phi);
1956 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
1957 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
1960 Double_t cellE = cells->GetAmplitude(i);
1961 cellIsolation += cellE;
1963 return cellIsolation;
1966 //________________________________________________________________________
1967 Double_t AliAnalysisTaskEMCALPi0PbPb::GetCellEnergy(const AliVCluster *cluster) const
1969 // Get maximum energy of attached cell.
1972 Int_t ncells = cluster->GetNCells();
1974 for (Int_t i=0; i<ncells; i++) {
1975 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1979 for (Int_t i=0; i<ncells; i++) {
1980 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1987 //________________________________________________________________________
1988 Double_t AliAnalysisTaskEMCALPi0PbPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1990 // Get maximum energy of attached cell.
1994 Int_t ncells = cluster->GetNCells();
1996 for (Int_t i=0; i<ncells; i++) {
1997 Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2000 id = cluster->GetCellAbsId(i);
2004 for (Int_t i=0; i<ncells; i++) {
2005 Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2008 id = cluster->GetCellAbsId(i);
2014 //________________________________________________________________________
2015 void AliAnalysisTaskEMCALPi0PbPb::GetSigma(const AliVCluster *c, Double_t& sigmaMax, Double_t &sigmaMin) const
2017 // Calculate the (E) weighted variance along the longer (eigen) axis.
2019 sigmaMax = 0; // cluster variance along its longer axis
2020 sigmaMin = 0; // cluster variance along its shorter axis
2021 Double_t Ec = c->E(); // cluster energy
2024 Double_t Xc = 0 ; // cluster first moment along X
2025 Double_t Yc = 0 ; // cluster first moment along Y
2026 Double_t Sxx = 0 ; // cluster second central moment along X (variance_X^2)
2027 Double_t Sxy = 0 ; // cluster second central moment along Y (variance_Y^2)
2028 Double_t Syy = 0 ; // cluster covariance^2
2030 AliVCaloCells *cells = fEsdCells;
2037 Int_t ncells = c->GetNCells();
2041 for(Int_t j=0; j<ncells; ++j) {
2042 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2043 Double_t cellen = cells->GetCellAmplitude(id);
2045 fGeom->GetGlobal(id,pos);
2046 Xc += cellen*pos.X();
2047 Yc += cellen*pos.Y();
2048 Sxx += cellen*pos.X()*pos.X();
2049 Syy += cellen*pos.Y()*pos.Y();
2050 Sxy += cellen*pos.X()*pos.Y();
2060 Sxx = TMath::Abs(Sxx);
2061 Syy = TMath::Abs(Syy);
2062 sigmaMax = (Sxx + Syy + TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2063 sigmaMax = TMath::Sqrt(TMath::Abs(sigmaMax));
2064 sigmaMin = TMath::Abs(Sxx + Syy - TMath::Sqrt(TMath::Abs((Sxx-Syy)*(Sxx-Syy)+4.0*Sxy*Sxy)))/2.0;
2065 sigmaMin = TMath::Sqrt(TMath::Abs(sigmaMin));
2068 //________________________________________________________________________
2069 Int_t AliAnalysisTaskEMCALPi0PbPb::GetNCells(const AliVCluster *c, Double_t emin) const
2071 // Calculate number of attached cells above emin.
2073 AliVCaloCells *cells = fEsdCells;
2080 Int_t ncells = c->GetNCells();
2081 for(Int_t j=0; j<ncells; ++j) {
2082 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2083 Double_t cellen = cells->GetCellAmplitude(id);
2090 //________________________________________________________________________
2091 Double_t AliAnalysisTaskEMCALPi0PbPb::GetTrackIsolation(Double_t cEta, Double_t cPhi, Double_t radius, Double_t pt) const
2093 // Compute isolation based on tracks.
2095 Double_t trkIsolation = 0;
2096 Double_t rad2 = radius*radius;
2097 Int_t ntrks = fSelPrimTracks->GetEntries();
2098 for(Int_t j = 0; j<ntrks; ++j) {
2099 AliVTrack *track = static_cast<AliVTrack*>(fSelTracks->At(j));
2104 Float_t eta = track->Eta();
2105 Float_t phi = track->Phi();
2106 Double_t phidiff = TVector2::Phi_mpi_pi(phi-cPhi);
2107 Double_t dist = (eta-cEta)*(eta-cEta)+phidiff*phidiff;
2110 trkIsolation += track->Pt();
2112 return trkIsolation;
2115 //________________________________________________________________________
2116 Bool_t AliAnalysisTaskEMCALPi0PbPb::IsShared(const AliVCluster *c) const
2118 // Returns if cluster shared across super modules.
2120 AliVCaloCells *cells = fEsdCells;
2127 Int_t ncells = c->GetNCells();
2128 for(Int_t j=0; j<ncells; ++j) {
2129 Int_t id = TMath::Abs(c->GetCellAbsId(j));
2130 Int_t got = id / (24*48);
2141 //________________________________________________________________________
2142 void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliVParticle *p, const TObjArray *arr, Int_t level) const
2144 // Print recursively daughter information.
2149 const AliAODMCParticle *amc = dynamic_cast<const AliAODMCParticle*>(p);
2152 for (Int_t i=0; i<level; ++i) printf(" ");
2155 Int_t n = amc->GetNDaughters();
2156 for (Int_t i=0; i<n; ++i) {
2157 Int_t d = amc->GetDaughter(i);
2158 const AliVParticle *dmc = static_cast<const AliVParticle*>(arr->At(d));
2159 PrintDaughters(dmc,arr,level+1);
2163 //________________________________________________________________________
2164 void AliAnalysisTaskEMCALPi0PbPb::PrintDaughters(const AliMCParticle *p, const AliMCEvent *arr, Int_t level) const
2166 // Print recursively daughter information.
2171 for (Int_t i=0; i<level; ++i) printf(" ");
2172 Int_t d1 = p->GetFirstDaughter();
2173 Int_t d2 = p->GetLastDaughter();
2174 printf("pid=%d: %.2f %.2f %.2f (%.2f %.2f %.2f); nd=%d,%d\n",
2175 p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2);
2180 for (Int_t i=d1;i<=d2;++i) {
2181 const AliMCParticle *dmc = static_cast<const AliMCParticle *>(arr->GetTrack(i));
2182 PrintDaughters(dmc,arr,level+1);
2186 //________________________________________________________________________
2187 void AliAnalysisTaskEMCALPi0PbPb::PrintTrackRefs(AliMCParticle *p) const
2189 // Print track ref array.
2194 Int_t n = p->GetNumberOfTrackReferences();
2195 for (Int_t i=0; i<n; ++i) {
2196 AliTrackReference *ref = p->GetTrackReference(i);
2199 ref->SetUserId(ref->DetectorId());
2204 //________________________________________________________________________
2205 void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliVParticle *p, Int_t index, const TObjArray *arr)
2207 // Process and create daughters.
2212 AliAODMCParticle *amc = dynamic_cast<AliAODMCParticle*>(p);
2218 Int_t nparts = arr->GetEntries();
2219 Int_t nents = fMcParts->GetEntries();
2221 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2222 newp->fPt = amc->Pt();
2223 newp->fEta = amc->Eta();
2224 newp->fPhi = amc->Phi();
2225 if (amc->Xv() != 0 || amc->Yv() != 0 || amc->Zv() != 0) {
2226 TVector3 vec(amc->Xv(),amc->Yv(),amc->Zv());
2227 newp->fVR = vec.Perp();
2228 newp->fVEta = vec.Eta();
2229 newp->fVPhi = vec.Phi();
2231 newp->fPid = amc->PdgCode();
2233 Int_t moi = amc->GetMother();
2234 if (moi>=0&&moi<nparts) {
2235 const AliAODMCParticle *mmc = static_cast<const AliAODMCParticle*>(arr->At(moi));
2236 moi = mmc->GetUniqueID();
2239 p->SetUniqueID(nents);
2241 // TODO: Determine which detector was hit
2244 Int_t n = amc->GetNDaughters();
2245 for (Int_t i=0; i<n; ++i) {
2246 Int_t d = amc->GetDaughter(i);
2247 if (d<=index || d>=nparts)
2249 AliVParticle *dmc = static_cast<AliVParticle*>(arr->At(d));
2250 ProcessDaughters(dmc,d,arr);
2254 //________________________________________________________________________
2255 void AliAnalysisTaskEMCALPi0PbPb::ProcessDaughters(AliMCParticle *p, Int_t index, const AliMCEvent *arr)
2257 // Process and create daughters.
2262 Int_t d1 = p->GetFirstDaughter();
2263 Int_t d2 = p->GetLastDaughter();
2265 printf("%d pid=%d: %.3f %.3f %.3f (%.2f %.2f %.2f); nd=%d,%d, mo=%d\n",
2266 index,p->PdgCode(),p->Px(),p->Py(),p->Pz(),p->Xv(),p->Yv(),p->Zv(),d1,d2, p->GetMother());
2269 Int_t nents = fMcParts->GetEntries();
2271 AliStaPart *newp = static_cast<AliStaPart*>(fMcParts->New(nents));
2272 newp->fPt = p->Pt();
2273 newp->fEta = p->Eta();
2274 newp->fPhi = p->Phi();
2275 if (p->Xv() != 0 || p->Yv() != 0 || p->Zv() != 0) {
2276 TVector3 vec(p->Xv(),p->Yv(),p->Zv());
2277 newp->fVR = vec.Perp();
2278 newp->fVEta = vec.Eta();
2279 newp->fVPhi = vec.Phi();
2281 newp->fPid = p->PdgCode();
2283 Int_t moi = p->GetMother();
2285 const AliMCParticle *mmc = static_cast<const AliMCParticle *>(arr->GetTrack(moi));
2286 moi = mmc->GetUniqueID();
2289 p->SetUniqueID(nents);
2291 Int_t nref = p->GetNumberOfTrackReferences();
2293 AliTrackReference *ref = p->GetTrackReference(nref-1);
2295 newp->fDet = ref->DetectorId();
2303 for (Int_t i=d1;i<=d2;++i) {
2304 AliMCParticle *dmc = static_cast<AliMCParticle *>(arr->GetTrack(i));
2307 ProcessDaughters(dmc,i,arr);