PWGGA/EMCALTasks/AliAnalysisTaskEMCALPi0PbPb.cxx

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