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