PWG4/CaloCalib/AliAnalysisTaskEMCALPi0PbPb.cxx

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