PWGHF/hfe/AliAnalysisTaskHFECal.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // Class for heavy-flavour electron  with EMCal triggered events
  17 // Author: Shingo Sakai
  18
  19
  20 #include "TChain.h"
  21 #include "TTree.h"
  22 #include "TH2F.h"
  23 #include "TF1.h"
  24 #include "TMath.h"
  25 #include "TCanvas.h"
  26 #include "THnSparse.h"
  27 #include "TLorentzVector.h"
  28 #include "TString.h"
  29 #include "TFile.h"
  30 #include "TGraphErrors.h"
  31
  32 #include "TDatabasePDG.h"
  33
  34 #include "AliAnalysisTask.h"
  35 #include "AliAnalysisManager.h"
  36
  37 #include "AliESDEvent.h"
  38 #include "AliESDHandler.h"
  39 #include "AliAODEvent.h"
  40 #include "AliAODHandler.h"
  41
  42 #include "AliMCEventHandler.h"
  43 #include "AliMCEvent.h"
  44 #include "AliMCParticle.h"
  45
  46 #include "AliAnalysisTaskHFECal.h"
  47 #include "TGeoGlobalMagField.h"
  48 #include "AliLog.h"
  49 #include "AliAnalysisTaskSE.h"
  50 #include "TRefArray.h"
  51 #include "TVector.h"
  52 #include "AliESDInputHandler.h"
  53 #include "AliESDpid.h"
  54 #include "AliESDtrackCuts.h"
  55 #include "AliPhysicsSelection.h"
  56 #include "AliESDCaloCluster.h"
  57 #include "AliAODCaloCluster.h"
  58 #include "AliEMCALRecoUtils.h"
  59 #include "AliEMCALGeometry.h"
  60 #include "AliGeomManager.h"
  61 #include "stdio.h"
  62 #include "TGeoManager.h"
  63 #include "iostream"
  64 #include "fstream"
  65
  66 #include "AliEMCALTrack.h"
  67 #include "AliMagF.h"
  68
  69 #include "AliKFParticle.h"
  70 #include "AliKFVertex.h"
  71
  72 #include "AliPID.h"
  73 #include "AliPIDResponse.h"
  74 #include "AliHFEcontainer.h"
  75 #include "AliHFEcuts.h"
  76 #include "AliHFEpid.h"
  77 #include "AliHFEpidBase.h"
  78 #include "AliHFEpidQAmanager.h"
  79 #include "AliHFEtools.h"
  80 #include "AliCFContainer.h"
  81 #include "AliCFManager.h"
  82
  83 #include "AliStack.h"
  84
  85 #include "AliCentrality.h"
  86
  87 using namespace std;
  88
  89 ClassImp(AliAnalysisTaskHFECal)
  90 //________________________________________________________________________
  91 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
  92   : AliAnalysisTaskSE(name)
  93   ,fESD(0)
  94   ,fMC(0)
  95   ,stack(0)
  96   ,fGeom(0)
  97   ,fOutputList(0)
  98   ,fqahist(0)
  99   ,fTrackCuts(0)
 100   ,fCuts(0)
 101   ,fIdentifiedAsOutInz(kFALSE)
 102   ,fPassTheEventCut(kFALSE)
 103   ,fRejectKinkMother(kFALSE)
 104   ,fmcData(kFALSE)
 105   ,fVz(0.0)
 106   ,fCFM(0)
 107   ,fPID(0)
 108   ,fPIDqa(0)
 109   ,fOpeningAngleCut(0.1)
 110   ,fMimpTassCut(0.5)
 111   ,fMimNsigassCut(-3)
 112   ,fInvmassCut(0)        // no mass
 113   ,fSetMassConstraint(kTRUE)
 114   ,fSetMassWidthCut(kFALSE)
 115   ,fSetMassNonlinear(kFALSE)
 116   ,fSetKFpart(kTRUE)
 117   ,fNoEvents(0)
 118   ,fEMCAccE(0)
 119   ,hEMCAccE(0)
 120   ,fTrkpt(0)
 121   ,fTrkEovPBef(0)
 122   ,fTrkEovPAft(0)
 123   ,fdEdxBef(0)
 124   ,fdEdxAft(0)
 125   ,fIncpT(0)
 126   ,fIncpTM20(0)
 127   ,fInvmassLS(0)
 128   ,fInvmassULS(0)
 129   ,fInvmassLSmc(0)
 130   ,fInvmassULSmc(0)
 131   ,fInvmassLSreco(0)
 132   ,fInvmassULSreco(0)
 133   ,fInvmassLSmc0(0)
 134   ,fInvmassLSmc1(0)
 135   ,fInvmassLSmc2(0)
 136   ,fInvmassLSmc3(0)
 137   ,fInvmassULSmc0(0)
 138   ,fInvmassULSmc1(0)
 139   ,fInvmassULSmc2(0)
 140   ,fInvmassULSmc3(0)
 141   ,fOpeningAngleLS(0)
 142   ,fOpeningAngleULS(0)
 143   ,fPhotoElecPt(0)
 144   ,fPhoElecPt(0)
 145   ,fPhoElecPtM20(0)
 146   ,fPhoElecPtM20Mass(0)
 147   ,fSameElecPt(0)
 148   ,fSameElecPtM20(0)
 149   ,fSameElecPtM20Mass(0)
 150   ,fSemiElecPtM20(0)
 151   ,fTrackPtBefTrkCuts(0)
 152   ,fTrackPtAftTrkCuts(0)
 153   ,fTPCnsigma(0)
 154   ,fCent(0)
 155   ,fEleInfo(0)
 156   ,fElenSigma(0)
 157   /*
 158   ,fClsEBftTrigCut(0)
 159   ,fClsEAftTrigCut(0)
 160   ,fClsEAftTrigCut1(0)
 161   ,fClsEAftTrigCut2(0)
 162   ,fClsEAftTrigCut3(0)
 163   ,fClsEAftTrigCut4(0)
 164   ,fClsETime(0)
 165   ,fClsETime1(0)
 166   ,fTrigTimes(0)
 167   ,fCellCheck(0)
 168   */
 169   ,fInputHFEMC(0)
 170   ,fInputAlle(0)
 171   ,fIncpTMChfe(0)
 172   ,fIncpTMChfeAll(0)
 173   ,fIncpTMCM20hfe(0)
 174   ,fIncpTMCM20hfeAll(0)
 175   ,fIncpTMCM20hfeCheck(0)
 176   ,fInputHFEMC_weight(0)
 177   ,fIncpTMCM20hfeCheck_weight(0)
 178   ,fIncpTMCpho(0)
 179   ,fIncpTMCM20pho(0)
 180   ,fPhoElecPtMC(0)
 181   ,fPhoElecPtMCM20(0)
 182   ,fPhoElecPtMCM20Mass(0)
 183   ,fSameElecPtMC(0)
 184   ,fSameElecPtMCM20(0)
 185   ,fSameElecPtMCM20Mass(0)
 186   ,fIncpTMCM20pho_pi0e(0)
 187   ,fPhoElecPtMCM20_pi0e(0)
 188   ,fSameElecPtMCM20_pi0e(0)
 189   ,fIncpTMCM20pho_eta(0)
 190   ,fPhoElecPtMCM20_eta(0)
 191   ,fSameElecPtMCM20_eta(0)
 192   ,fIncpTMCpho_pi0e_TPC(0)
 193   ,fPhoElecPtMC_pi0e_TPC(0)
 194   ,fSameElecPtMC_pi0e_TPC(0)
 195   ,fIncpTMCpho_eta_TPC(0)
 196   ,fPhoElecPtMC_eta_TPC(0)
 197   ,fSameElecPtMC_eta_TPC(0)
 198   ,CheckNclust(0)
 199   ,CheckNits(0)
 200   ,CheckDCA(0)
 201   ,Hpi0pTcheck(0)
 202   ,HETApTcheck(0)
 203   ,HphopTcheck(0)
 204   ,HDpTcheck(0)
 205   ,HBpTcheck(0)
 206   ,fpTCheck(0)
 207   ,fMomDtoE(0)
 208   ,fLabelCheck(0)
 209   ,fgeoFake(0)
 210   ,fFakeTrk0(0)
 211   ,fFakeTrk1(0)
 212   ,ftimingEle(0)
 213   ,fIncMaxE(0)
 214   ,fIncReco(0)
 215   ,fPhoReco(0)
 216   ,fSamReco(0)
 217   ,fIncRecoMaxE(0)
 218   ,fPhoRecoMaxE(0)
 219   ,fSamRecoMaxE(0)
 220   ,fPhoVertexReco_TPC(0)
 221   ,fPhoVertexReco_TPC_Invmass(0)
 222   ,fPhoVertexReco_EMCal(0)
 223   ,fPhoVertexReco_Invmass(0)
 224   ,fPhoVertexReco_step0(0)
 225   ,fPhoVertexReco_step1(0)
 226   ,fMatchV0_0(0)
 227   ,fMatchV0_1(0)
 228   ,fMatchMC_0(0)
 229   ,fMatchMC_1(0)
 230   ,fpair(0)
 231   ,fFakeRejection0(0)
 232   ,fFakeRejection1(0)
 233   ,fFakeRejection2(0)
 234   ,EopFake(0)
 235   ,EopTrue(0)
 236   ,MatchFake(0)
 237   ,MatchTrue(0)
 238   //,fnSigEtaCorr(NULL)
 239 {
 240   //Named constructor
 241
 242   fPID = new AliHFEpid("hfePid");
 243   fTrackCuts = new AliESDtrackCuts();
 244
 245   for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0;
 246
 247   // Define input and output slots here
 248   // Input slot #0 works with a TChain
 249   DefineInput(0, TChain::Class());
 250   // Output slot #0 id reserved by the base class for AOD
 251   // Output slot #1 writes into a TH1 container
 252   // DefineOutput(1, TH1I::Class());
 253   DefineOutput(1, TList::Class());
 254   //  DefineOutput(3, TTree::Class());
 255 }
 256
 257 //________________________________________________________________________
 258 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
 259   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
 260   ,fESD(0)
 261   ,fMC(0)
 262   ,stack(0)
 263   ,fGeom(0)
 264   ,fOutputList(0)
 265   ,fqahist(0)
 266   ,fTrackCuts(0)
 267   ,fCuts(0)
 268   ,fIdentifiedAsOutInz(kFALSE)
 269   ,fPassTheEventCut(kFALSE)
 270   ,fRejectKinkMother(kFALSE)
 271   ,fmcData(kFALSE)
 272   ,fVz(0.0)
 273   ,fCFM(0)
 274   ,fPID(0)
 275   ,fPIDqa(0)
 276   ,fOpeningAngleCut(0.1)
 277   ,fMimpTassCut(0.5)
 278   ,fMimNsigassCut(-3)
 279   ,fInvmassCut(0)        // no mass
 280   ,fSetMassConstraint(kTRUE)
 281   ,fSetMassWidthCut(kFALSE)
 282   ,fSetMassNonlinear(kFALSE)
 283   ,fSetKFpart(kTRUE)
 284   ,fNoEvents(0)
 285   ,fEMCAccE(0)
 286   ,hEMCAccE(0)
 287   ,fTrkpt(0)
 288   ,fTrkEovPBef(0)
 289   ,fTrkEovPAft(0)
 290   ,fdEdxBef(0)
 291   ,fdEdxAft(0)
 292   ,fIncpT(0)
 293   ,fIncpTM20(0)
 294   ,fInvmassLS(0)
 295   ,fInvmassULS(0)
 296   ,fInvmassLSmc(0)
 297   ,fInvmassULSmc(0)
 298   ,fInvmassLSreco(0)
 299   ,fInvmassULSreco(0)
 300   ,fInvmassLSmc0(0)
 301   ,fInvmassLSmc1(0)
 302   ,fInvmassLSmc2(0)
 303   ,fInvmassLSmc3(0)
 304   ,fInvmassULSmc0(0)
 305   ,fInvmassULSmc1(0)
 306   ,fInvmassULSmc2(0)
 307   ,fInvmassULSmc3(0)
 308   ,fOpeningAngleLS(0)
 309   ,fOpeningAngleULS(0)
 310   ,fPhotoElecPt(0)
 311   ,fPhoElecPt(0)
 312   ,fPhoElecPtM20(0)
 313   ,fPhoElecPtM20Mass(0)
 314   ,fSameElecPt(0)
 315   ,fSameElecPtM20(0)
 316   ,fSameElecPtM20Mass(0)
 317   ,fSemiElecPtM20(0)
 318   ,fTrackPtBefTrkCuts(0)
 319   ,fTrackPtAftTrkCuts(0)
 320   ,fTPCnsigma(0)
 321   ,fCent(0)
 322   ,fEleInfo(0)
 323   ,fElenSigma(0)
 324   /*
 325   ,fClsEBftTrigCut(0)
 326   ,fClsEAftTrigCut(0)
 327   ,fClsEAftTrigCut1(0)
 328   ,fClsEAftTrigCut2(0)
 329   ,fClsEAftTrigCut3(0)
 330   ,fClsEAftTrigCut4(0)
 331   ,fClsETime(0)
 332   ,fClsETime1(0)
 333   ,fTrigTimes(0)
 334   ,fCellCheck(0)
 335   */
 336   ,fInputHFEMC(0)
 337   ,fInputAlle(0)
 338   ,fIncpTMChfe(0)
 339   ,fIncpTMChfeAll(0)
 340   ,fIncpTMCM20hfe(0)
 341   ,fIncpTMCM20hfeAll(0)
 342   ,fIncpTMCM20hfeCheck(0)
 343   ,fInputHFEMC_weight(0)
 344   ,fIncpTMCM20hfeCheck_weight(0)
 345   ,fIncpTMCpho(0)
 346   ,fIncpTMCM20pho(0)
 347   ,fPhoElecPtMC(0)
 348   ,fPhoElecPtMCM20(0)
 349   ,fPhoElecPtMCM20Mass(0)
 350   ,fSameElecPtMC(0)
 351   ,fSameElecPtMCM20(0)
 352   ,fSameElecPtMCM20Mass(0)
 353   ,fIncpTMCM20pho_pi0e(0)
 354   ,fPhoElecPtMCM20_pi0e(0)
 355   ,fSameElecPtMCM20_pi0e(0)
 356   ,fIncpTMCM20pho_eta(0)
 357   ,fPhoElecPtMCM20_eta(0)
 358   ,fSameElecPtMCM20_eta(0)
 359   ,fIncpTMCpho_pi0e_TPC(0)
 360   ,fPhoElecPtMC_pi0e_TPC(0)
 361   ,fSameElecPtMC_pi0e_TPC(0)
 362   ,fIncpTMCpho_eta_TPC(0)
 363   ,fPhoElecPtMC_eta_TPC(0)
 364   ,fSameElecPtMC_eta_TPC(0)
 365   ,CheckNclust(0)
 366   ,CheckNits(0)
 367   ,CheckDCA(0)
 368   ,Hpi0pTcheck(0)
 369   ,HETApTcheck(0)
 370   ,HphopTcheck(0)
 371   ,HDpTcheck(0)
 372   ,HBpTcheck(0)
 373   ,fpTCheck(0)
 374   ,fMomDtoE(0)
 375   ,fLabelCheck(0)
 376   ,fgeoFake(0)
 377   ,fFakeTrk0(0)
 378   ,fFakeTrk1(0)
 379   ,ftimingEle(0)
 380   ,fIncMaxE(0)
 381   ,fIncReco(0)
 382   ,fPhoReco(0)
 383   ,fSamReco(0)
 384   ,fIncRecoMaxE(0)
 385   ,fPhoRecoMaxE(0)
 386   ,fSamRecoMaxE(0)
 387   ,fPhoVertexReco_TPC(0)
 388   ,fPhoVertexReco_TPC_Invmass(0)
 389   ,fPhoVertexReco_EMCal(0)
 390   ,fPhoVertexReco_Invmass(0)
 391   ,fPhoVertexReco_step0(0)
 392   ,fPhoVertexReco_step1(0)
 393   ,fMatchV0_0(0)
 394   ,fMatchV0_1(0)
 395   ,fMatchMC_0(0)
 396   ,fMatchMC_1(0)
 397   ,fpair(0)
 398   ,fFakeRejection0(0)
 399   ,fFakeRejection1(0)
 400   ,fFakeRejection2(0)
 401   ,EopFake(0)
 402   ,EopTrue(0)
 403   ,MatchFake(0)
 404   ,MatchTrue(0)
 405   //,fnSigEtaCorr(NULL)
 406 {
 407         //Default constructor
 408         fPID = new AliHFEpid("hfePid");
 409
 410         fTrackCuts = new AliESDtrackCuts();
 411
 412        for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0;
 413
 414         // Constructor
 415         // Define input and output slots here
 416         // Input slot #0 works with a TChain
 417         DefineInput(0, TChain::Class());
 418         // Output slot #0 id reserved by the base class for AOD
 419         // Output slot #1 writes into a TH1 container
 420         // DefineOutput(1, TH1I::Class());
 421         DefineOutput(1, TList::Class());
 422         //DefineOutput(3, TTree::Class());
 423 }
 424 //_________________________________________
 425
 426 AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
 427 {
 428   //Destructor
 429
 430   delete fOutputList;
 431   delete fGeom;
 432   delete fPID;
 433   delete fCuts;
 434   delete fCFM;
 435   delete fPIDqa;
 436   delete fTrackCuts;
 437 }
 438 //_________________________________________
 439
 440 void AliAnalysisTaskHFECal::UserExec(Option_t*)
 441 {
 442   //Main loop
 443   //Called for each event
 444
 445   // create pointer to event
 446   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 447   if (!fESD) {
 448     printf("ERROR: fESD not available\n");
 449     return;
 450   }
 451
 452   if(!fCuts){
 453     AliError("HFE cuts not available");
 454     return;
 455   }
 456
 457   if(!fPID->IsInitialized()){
 458     // Initialize PID with the given run number
 459     AliWarning("PID not initialised, get from Run no");
 460     fPID->InitializePID(fESD->GetRunNumber());
 461   }
 462
 463   if(fmcData)fMC = MCEvent();
 464   //AliStack* stack = NULL;
 465   if(fmcData && fMC)stack = fMC->Stack();
 466
 467   Float_t cent = -1.;
 468   AliCentrality *centrality = fESD->GetCentrality();
 469   cent = centrality->GetCentralityPercentile("V0M");
 470
 471   //---- fill MC track info
 472   if(fmcData && fMC)
 473     {
 474     Int_t nParticles = stack->GetNtrack();
 475     for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
 476       TParticle* particle = stack->Particle(iParticle);
 477       int fPDG = particle->GetPdgCode();
 478       double mcZvertex = fMC->GetPrimaryVertex()->GetZ();
 479       double pTMC = particle->Pt();
 480       double proR = particle->R();
 481       double etaMC = particle->Eta();
 482       if(fabs(etaMC)>0.6)continue;
 483       Bool_t mcInDtoE= kFALSE;
 484       Bool_t mcInBtoE= kFALSE;
 485
 486       Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
 487       //if(!MChijing)printf("not MC hijing");
 488       int iHijing = 1;
 489       if(!MChijing)iHijing = 0;
 490       double mcphoinfo[5];
 491       mcphoinfo[0] = cent;
 492       mcphoinfo[1] = pTMC;
 493       mcphoinfo[2] = iHijing;
 494       if(fPDG==111)Hpi0pTcheck->Fill(mcphoinfo);
 495       if(fPDG==221)HETApTcheck->Fill(mcphoinfo);
 496       if(fabs(fPDG)==411 || fabs(fPDG)==413 || fabs(fPDG)==421 || fabs(fPDG)==423 || fabs(fPDG)==431)HDpTcheck->Fill(pTMC,iHijing);
 497       if(fabs(fPDG)==511 || fabs(fPDG)==513 || fabs(fPDG)==521 || fabs(fPDG)==523 || fabs(fPDG)==531)HBpTcheck->Fill(pTMC,iHijing);
 498
 499       if(particle->GetFirstMother()>-1 && fPDG==22)
 500         {
 501          int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
 502          if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g
 503         }
 504
 505       if(particle->GetFirstMother()>-1 && fabs(fPDG)==11)
 506         {
 507             int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
 508             double pTMCparent = stack->Particle(particle->GetFirstMother())->Pt();
 509             if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE;
 510             if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE;
 511             if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0)
 512                {
 513                 fInputHFEMC->Fill(cent,pTMC);
 514                 double mcinfo[5];
 515                 mcinfo[0] = cent;
 516                 mcinfo[1] = pTMC;
 517                 mcinfo[2] = 0.0;
 518                 mcinfo[3] = iHijing;
 519                 mcinfo[4] = pTMCparent;
 520                 fInputHFEMC_weight->Fill(mcinfo);
 521                }
 522          }
 523
 524
 525          if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC);
 526
 527       }
 528     }
 529
 530   fNoEvents->Fill(0);
 531
 532   Int_t fNOtrks =  fESD->GetNumberOfTracks();
 533   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 534
 535   Double_t pVtxZ = -999;
 536   pVtxZ = pVtx->GetZ();
 537
 538   if(TMath::Abs(pVtxZ)>10) return;
 539   fNoEvents->Fill(1);
 540
 541   if(fNOtrks<1) return;
 542   fNoEvents->Fill(2);
 543
 544   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 545   if(!pidResponse){
 546     AliDebug(1, "Using default PID Response");
 547     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 548   }
 549
 550   fPID->SetPIDResponse(pidResponse);
 551
 552   fCFM->SetRecEventInfo(fESD);
 553
 554   fCent->Fill(cent);
 555
 556   //if(cent>90.) return;
 557
 558  // Calorimeter info.
 559
 560    FindTriggerClusters();
 561
 562   // make EMCAL array
 563   double maxE = 0.0;
 564   for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
 565      {
 566       AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
 567       if(clust && clust->IsEMCAL())
 568         {
 569          double clustE = clust->E();
 570          float  emcx[3]; // cluster pos
 571          clust->GetPosition(emcx);
 572          TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
 573          double emcphi = clustpos.Phi();
 574          double emceta = clustpos.Eta();
 575          double calInfo[5];
 576          calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2();
 577          //fEMCAccE->Fill(calInfo);
 578          hEMCAccE->Fill(cent,clustE);
 579          if(clustE>maxE)maxE = clustE;
 580         }
 581    }
 582
 583   // Track loop
 584   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 585     AliESDtrack* track = fESD->GetTrack(iTracks);
 586     if (!track) {
 587       printf("ERROR: Could not receive track %d\n", iTracks);
 588       continue;
 589     }
 590
 591     //--- Get MC informtion
 592
 593     int parentlabel = 99999;
 594     int parentPID = 99999;
 595     int grand_parentlabel = 99999;
 596     int grand_parentPID = 99999;
 597     Bool_t mcPho = kFALSE;
 598     Bool_t mcDtoE= kFALSE;
 599     Bool_t mcBtoE= kFALSE;
 600     Bool_t mcOrgPi0 = kFALSE;
 601     Bool_t mcOrgEta = kFALSE;
 602     double mcele = -1.;
 603     double mcpT = 0.0;
 604     double mcMompT = 0.0;
 605     //double mcGrandMompT = 0.0;
 606     double mcWeight = -10.0;
 607     double conv_proR = -1.0;
 608
 609     int iHijing = 1;
 610     int mcLabel = -1;
 611
 612     if(fmcData && fMC && stack)
 613       {
 614        Int_t label = TMath::Abs(track->GetLabel());
 615        //mcLabel = track->GetLabel();
 616        mcLabel = abs(track->GetLabel()); // check for conv. issue
 617
 618        if(mcLabel>-1)
 619        {
 620
 621                Bool_t MChijing = fMC->IsFromBGEvent(label);
 622                if(!MChijing)iHijing = 0;
 623
 624                TParticle* particle = stack->Particle(label);
 625                int mcpid = particle->GetPdgCode();
 626                mcpT = particle->Pt();
 627                conv_proR = particle->R();
 628                //printf("MCpid = %d",mcpid);
 629                if(particle->GetFirstMother()>-1)
 630                {
 631                        //int parentlabel = particle->GetFirstMother();
 632                        //int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode();
 633                        //mcMompT = stack->Particle(particle->GetFirstMother())->Pt();
 634                        FindMother(particle, parentlabel, parentPID);
 635                        mcMompT = stack->Particle(parentlabel)->Pt();
 636                        if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE;
 637                        if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE;
 638                        if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE;
 639
 640                        // make D->e pT correlation
 641                        if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT);
 642
 643                        //cout << "check PID = " << parentPID << endl;
 644                        //cout << "check pho = " << mcPho << endl;
 645                        //cout << "check D or B = " << mcDtoE << endl;
 646                        // pi->e (Dalitz)
 647                        if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0)
 648                        {
 649                                //cout << "find pi0->e " <<  endl;
 650                                mcOrgPi0 = kTRUE;
 651                                mcWeight = GetMCweight(mcMompT);
 652                        }
 653                        // eta->e (Dalitz)
 654                        if(parentPID==221 && fabs(mcpid)==11 && mcMompT>0.0)
 655                        {
 656                                //cout << "find Eta->e " <<  endl;
 657                                mcOrgEta = kTRUE;
 658                                mcWeight = GetMCweightEta(mcMompT);
 659                        }
 660
 661                        // access grand parent
 662                        TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer
 663                        //if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e
 664                        if(particle_parent->GetFirstMother()>-1 && (parentPID==22 || parentPID==111) && fabs(mcpid)==11) // get grand parent g->e
 665                        {
 666                                //int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode();
 667                                //double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt();
 668                                FindMother(particle_parent, grand_parentlabel, grand_parentPID);
 669                                double mcGrandpT = stack->Particle(grand_parentlabel)->Pt();
 670                                if(grand_parentPID==111 && mcGrandpT>0.0)
 671                                {
 672                                        // check eta->pi0 decay !
 673                                        int grand2_parentlabel = 99999; int grand2_parentPID = 99999;
 674                                        TParticle* particle_grand = stack->Particle(grand_parentlabel); // get parent pointer
 675                                        FindMother(particle_grand, grand2_parentlabel, grand2_parentPID);
 676                                        if(grand2_parentPID==221)
 677                                        {
 678                                                //cout << "find Eta->e " <<  endl;
 679                                                double mcGrandpT2 = stack->Particle(grand2_parentlabel)->Pt();
 680                                                mcOrgEta = kTRUE;
 681                                                mcWeight = GetMCweight(mcGrandpT2);
 682                                                mcMompT = mcGrandpT2;
 683                                        }
 684                                        else
 685                                        {
 686                                                //cout << "find pi0->e " <<  endl;
 687                                                mcOrgPi0 = kTRUE;
 688                                                mcWeight = GetMCweight(mcGrandpT);
 689                                                mcMompT = mcGrandpT;
 690                                        }
 691                                }
 692
 693                                if(grand_parentPID==221 && mcGrandpT>0.0)
 694                                {
 695                                        //cout << "find Eta->e " <<  endl;
 696                                        mcOrgEta = kTRUE;
 697                                        mcOrgPi0 = kFALSE;
 698                                        mcWeight = GetMCweightEta(mcGrandpT);
 699                                        mcMompT = mcGrandpT;
 700                                }
 701                        }
 702                }
 703
 704                //cout << "===================="<<endl;
 705                //cout << "mcDtoE : " << mcDtoE << endl;
 706                //cout << "mcBtoE : " << mcBtoE << endl;
 707                //cout << "mcPho : " << mcPho << endl;
 708
 709                if(fabs(mcpid)==11)mcele= 0.;
 710                //cout << "check e: " << mcele << endl;
 711                if(fabs(mcpid)==11 && mcDtoE)mcele= 1.;
 712                //cout << "check D->e: " << mcele << endl;
 713                if(fabs(mcpid)==11 && mcBtoE)mcele= 2.;
 714                //cout << "check B->e: " << mcele << endl;
 715                if(fabs(mcpid)==11 && mcPho)mcele= 3.;
 716                //cout << "check Pho->e: " << mcele << endl;
 717
 718                //cout << "check PID " << endl;
 719                if(fabs(mcpid)!=11)
 720                  {
 721                   //cout << "!= 11" << endl;
 722                   //cout << mcpid << endl;
 723                  }
 724                if(mcele==-1)
 725                  {
 726                   //cout << "mcele==-1" << endl;
 727                   //cout << mcele << endl;
 728                   //cout << mcpid << endl;
 729                  }
 730
 731        } // end of mcLabel>-1
 732
 733       } // <------ end of MC info.
 734
 735     //cout << "Pi0 = " << mcOrgPi0 << " ; Eta = " << mcOrgEta << endl;
 736     //printf("weight = %f\n",mcWeight);
 737
 738     if(TMath::Abs(track->Eta())>0.6) continue;
 739     //if(TMath::Abs(track->Pt()<2.5)) continue;
 740     if(TMath::Abs(track->Pt()<0.1)) continue;
 741
 742     int nITS = track->GetNcls(0);
 743
 744     fTrackPtBefTrkCuts->Fill(track->Pt());
 745
 746     /*
 747     UChar_t itsPixel = track->GetITSClusterMap();
 748     cout << "nITS = " << nITS << endl;
 749     if(itsPixel & BIT(0))cout << "1st layer hit" << endl;
 750     if(itsPixel & BIT(1))cout << "2nd layer hit" << endl;
 751     if(itsPixel & BIT(2))cout << "3rd layer hit" << endl;
 752     if(itsPixel & BIT(3))cout << "4th layer hit" << endl;
 753     if(itsPixel & BIT(4))cout << "5th layer hit" << endl;
 754     */
 755
 756     // RecKine: ITSTPC cuts
 757     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 758
 759     //RecKink
 760     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 761       if(track->GetKinkIndex(0) != 0) continue;
 762     }
 763
 764     // RecPrim
 765     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 766
 767     // HFEcuts: ITS layers cuts
 768     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 769
 770     // HFE cuts: TPC PID cleanup
 771     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 772
 773     if(mcPho && iHijing==0)fPhoVertexReco_step0->Fill(track->Pt(),conv_proR); // check MC vertex
 774     if(mcPho && iHijing==1)fPhoVertexReco_step1->Fill(track->Pt(),conv_proR); // check MC vertex
 775
 776     int nTPCcl = track->GetTPCNcls();
 777     //int nTPCclF = track->GetTPCNclsF(); // warnings
 778     //int nITS = track->GetNcls(0);
 779
 780
 781     fTrackPtAftTrkCuts->Fill(track->Pt());
 782
 783     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
 784     pt = track->Pt();
 785     if(pt<2.5)continue;
 786     if(pt<0.1)continue;
 787
 788     //Int_t charge = track->Charge();
 789     fTrkpt->Fill(pt);
 790     mom = track->P();
 791     phi = track->Phi();
 792     eta = track->Eta();
 793     float dca_xy;
 794     float dca_z;
 795     track->GetImpactParameters(dca_xy,dca_z);
 796
 797     dEdx = track->GetTPCsignal();
 798     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 799
 800     //cout << "nSigma correctoon-----" << endl;
 801     //cout << "org = " << fTPCnSigma << endl;
 802     /*
 803     if(!fmcData) // nsigma eta correction
 804        {
 805         double nSigexpCorr = NsigmaCorrection(eta,cent);
 806         fTPCnSigma -= nSigexpCorr;
 807        }
 808     */
 809     //cout << "correction = " << fTPCnSigma << endl;
 810
 811     //--- track cluster match
 812
 813         double ncells = -1.0;
 814         double m20 = -1.0;
 815         double m02 = -1.0;
 816         double disp = -1.0;
 817         double rmatch = -1.0;
 818         double nmatch = -1.0;
 819         //double oppstatus = 0.0;
 820         double emctof = 0.0;
 821         Bool_t MaxEmatch = kFALSE;
 822
 823     Int_t clsId = track->GetEMCALcluster();
 824     if (clsId>0){
 825       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
 826       if(clust && clust->IsEMCAL()){
 827
 828                 double clustE = clust->E();
 829                 if(clustE==maxE)MaxEmatch = kTRUE;
 830                 eop = clustE/fabs(mom);
 831
 832                 //double clustT = clust->GetTOF();
 833                 ncells = clust->GetNCells();
 834                 m02 = clust->GetM02();
 835                 m20 = clust->GetM20();
 836                 disp = clust->GetDispersion();
 837                 double delphi = clust->GetTrackDx();
 838                 double deleta = clust->GetTrackDz();
 839                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
 840                 nmatch = clust->GetNTracksMatched();
 841                 emctof = clust->GetTOF();
 842                 //cout << "emctof = " << emctof << endl;
 843                 cout << "eop org = "<< eop << endl;
 844                 double eoporg =  eop;
 845                 if(fmcData)
 846                   {
 847                    double mceopcorr = MCEopMeanCorrection(pt,cent);
 848                    eop += mceopcorr;
 849                   }
 850                 cout << "eop corr = " << eop << endl;
 851
 852                   double valdedx[16];
 853                   valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
 854                   valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = dca_xy,  valdedx[8] = dca_z; valdedx[9] = m20; valdedx[10] = mcpT;
 855                   valdedx[11] = cent; valdedx[12] = dEdx; valdedx[13] = eoporg; valdedx[14] = nTPCcl;
 856                   valdedx[15] = mcele;
 857                   fEleInfo->Fill(valdedx);
 858
 859       }
 860     }
 861
 862     //Get Cal info PID response
 863     /*
 864     double eop2;
 865     double ss[4];
 866     Double_t nSigmaEop = fPID->GetPIDResponse()->NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss);
 867     if(fTPCnSigma>-1.5 && fTPCnSigma<3.0 && nITS>2.5 && nTPCcl>100)
 868       {
 869        double valEop[3];
 870        valEop[0] = cent;
 871        valEop[1] = pt;
 872        valEop[2] = nSigmaEop;
 873        fElenSigma->Fill(valEop);
 874       }
 875     */
 876    // --- tarck cut & e ID
 877
 878     if(nITS<2.5)continue;
 879     if(nTPCcl<100)continue;
 880
 881
 882     CheckNclust->Fill(nTPCcl);
 883     CheckNits->Fill(nITS);
 884     CheckDCA->Fill(dca_xy,dca_z);
 885     // check production vertex of photons
 886
 887
 888     fdEdxBef->Fill(mom,fTPCnSigma);
 889     fTPCnsigma->Fill(mom,fTPCnSigma);
 890     if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
 891
 892     //--- track accepted by HFE
 893
 894     Int_t pidpassed = 1;
 895     AliHFEpidObject hfetrack;
 896     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 897     hfetrack.SetRecTrack(track);
 898     double binct = 10.5;
 899     if((0.0< cent) && (cent<5.0)) binct = 0.5;
 900     if((5.0< cent) && (cent<10.0)) binct = 1.5;
 901     if((10.0< cent) && (cent<20.0)) binct = 2.5;
 902     if((20.0< cent) && (cent<30.0)) binct = 3.5;
 903     if((30.0< cent) && (cent<40.0)) binct = 4.5;
 904     if((40.0< cent) && (cent<50.0)) binct = 5.5;
 905     if((50.0< cent) && (cent<60.0)) binct = 6.5;
 906     if((60.0< cent) && (cent<70.0)) binct = 7.5;
 907     if((70.0< cent) && (cent<80.0)) binct = 8.5;
 908     if((80.0< cent) && (cent<90.0)) binct = 9.5;
 909     if((90.0< cent) && (cent<100.0)) binct = 10.5;
 910
 911     hfetrack.SetCentrality((int)binct); //added
 912     hfetrack.SetPbPb();
 913     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 914
 915     if(pidpassed==0) continue; // nSigma rejection
 916
 917     //--- photonic ID
 918
 919     Bool_t fFlagPhotonicTPC = kFALSE;
 920     Bool_t fFlagConvinatTPC = kFALSE;
 921     SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicTPC,fFlagConvinatTPC,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta,0);
 922
 923     //--- check reco eff. with TPC
 924     double phoval[5];
 925     phoval[0] = cent;
 926     phoval[1] = pt;
 927     phoval[2] = fTPCnSigma;
 928     phoval[3] = iHijing;
 929     phoval[4] = mcMompT;
 930
 931     if((fTPCnSigma >= -5.0 && fTPCnSigma <= 5) && (mcOrgPi0 || mcOrgEta))
 932       {
 933         if(iHijing==1)mcWeight = 1.0;
 934         if(mcOrgPi0)
 935           {
 936            fIncpTMCpho_pi0e_TPC->Fill(phoval,mcWeight);
 937            if(fFlagPhotonicTPC) fPhoElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
 938            if(fFlagConvinatTPC) fSameElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
 939            fPhoVertexReco_TPC->Fill(track->Pt(),conv_proR,mcWeight); // check MC vertex
 940            if(fFlagPhotonicTPC)fPhoVertexReco_TPC_Invmass->Fill(track->Pt(),conv_proR,mcWeight); // check MC vertex
 941           }
 942         if(mcOrgEta)
 943           {
 944            fIncpTMCpho_eta_TPC->Fill(phoval,mcWeight);
 945            if(fFlagPhotonicTPC) fPhoElecPtMC_eta_TPC->Fill(phoval,mcWeight);
 946            if(fFlagConvinatTPC) fSameElecPtMC_eta_TPC->Fill(phoval,mcWeight);
 947           }
 948       }
 949
 950     //--- matching check
 951     double emcphimim = 1.396;
 952     double emcphimax = 3.14;
 953     if(phi>emcphimim && phi<emcphimax)
 954       {
 955        if(fFlagPhotonicTPC)fMatchV0_0->Fill(pt);  // data
 956        if(mcele>2.1)fMatchMC_0->Fill(pt,mcWeight); // MC
 957
 958        if(eop>=0.0) // have a match
 959          {
 960           if(fFlagPhotonicTPC)fMatchV0_1->Fill(pt);
 961           if(mcele>2.1)fMatchMC_1->Fill(pt,mcWeight);
 962          }
 963       }
 964
 965      // check fake rejection
 966     if(mcOrgPi0 || mcOrgEta)
 967       {
 968        double phiacc0 = 80.0/180.0*acos(-1);
 969        double phiacc1 = 180.0/180.0*acos(-1);
 970        int TrStat = 0;
 971        if(phi>phiacc0 && phi<phiacc1)
 972          {
 973           if(track->GetLabel()>0)TrStat = 1;
 974           fFakeRejection0->Fill(TrStat,pt,mcWeight);
 975           if(eop>-1.0)fFakeRejection1->Fill(TrStat,pt,mcWeight); // have match
 976           if(eop>0.9 && eop<1.3)fFakeRejection2->Fill(TrStat,pt,mcWeight); // have PID
 977
 978           if(TrStat==0)
 979             {
 980              EopFake->Fill(pt,eop);
 981              MatchFake->Fill(pt,rmatch);
 982             }
 983           else
 984             {
 985              EopTrue->Fill(pt,eop);
 986              MatchTrue->Fill(pt,rmatch);
 987             }
 988
 989         }
 990       }
 991
 992     //+++++++  E/p cut ++++++++++++++++
 993
 994     if(eop<0.9 || eop>1.3)continue;
 995
 996     Bool_t fFlagPhotonicElec = kFALSE;
 997     Bool_t fFlagConvinatElec = kFALSE;
 998     SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta,1);
 999
1000     fTrkEovPAft->Fill(pt,eop);
1001     fdEdxAft->Fill(mom,fTPCnSigma);
1002
1003     // Fill real data
1004     fIncpT->Fill(cent,pt);
1005     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
1006     if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
1007
1008     if(m20>0.0 && m20<0.3)
1009       {
1010        fIncpTM20->Fill(cent,pt);
1011        ftimingEle->Fill(pt,emctof);
1012        if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
1013        if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
1014        if(!fFlagPhotonicElec) fSemiElecPtM20->Fill(cent,pt);
1015      }
1016
1017
1018     // MC
1019     // check label for electron candidiates
1020     int idlabel = 1;
1021     if(mcLabel==0)idlabel = 0;
1022     fLabelCheck->Fill(pt,idlabel);
1023     if(mcLabel==0)fgeoFake->Fill(phi,eta);
1024
1025     if(mcLabel<0 && m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)
1026        {
1027         fFakeTrk0->Fill(cent,pt);
1028        }
1029
1030     if(mcele>-1) // select MC electrons
1031       {
1032
1033           fIncpTMChfeAll->Fill(cent,pt);
1034           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
1035           if(m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)fFakeTrk1->Fill(cent,pt);
1036
1037        if(mcBtoE || mcDtoE) // select B->e & D->e
1038          {
1039           fIncpTMChfe->Fill(cent,pt);
1040           if(m20>0.0 && m20<0.3)
1041             {
1042                 //cout << "MC label = " << mcLabel << endl;
1043                 fIncpTMCM20hfe->Fill(cent,pt);
1044                 fIncpTMCM20hfeCheck->Fill(cent,mcpT);
1045                 fIncpTMCM20hfeCheck_weight->Fill(phoval);
1046             }
1047          }
1048
1049        if(mcPho) // select photonic electrons
1050         {
1051
1052          fIncpTMCpho->Fill(phoval);
1053          if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
1054          if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
1055
1056          if(m20>0.0 && m20<0.3)
1057            {
1058             fIncpTMCM20pho->Fill(phoval);
1059             if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
1060             if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
1061             // pi0->g->e
1062             if(mcWeight>-1)
1063               {
1064                if(iHijing==1)mcWeight = 1.0;
1065                if(mcOrgPi0)
1066                  {
1067                   fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
1068                   if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
1069                   if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
1070
1071                   // check production vertex
1072                   //fPhoVertexReco_EMCal->Fill(track->Pt(),conv_proR);
1073                   //if(fFlagPhotonicElec)fPhoVertexReco_Invmass->Fill(track->Pt(),conv_proR);
1074                   fPhoVertexReco_EMCal->Fill(track->Pt(),conv_proR,mcWeight);
1075                   if(fFlagPhotonicElec)fPhoVertexReco_Invmass->Fill(track->Pt(),conv_proR,mcWeight);
1076
1077                   }
1078                // --- eta
1079                if(mcOrgEta)
1080                  {
1081                   fIncpTMCM20pho_eta->Fill(phoval,mcWeight);
1082                   if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
1083                   if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
1084                  }
1085                 // check production vertex
1086               }
1087            }
1088         }
1089       }
1090  }
1091  PostData(1, fOutputList);
1092 }
1093 //_________________________________________
1094 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
1095 {
1096   //--- Check MC
1097
1098   //Bool_t mcData = kFALSE;
1099   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
1100     {
1101      fmcData = kTRUE;
1102      printf("+++++ MC Data available");
1103     }
1104   if(fmcData)
1105     {
1106      printf("++++++++= MC analysis \n");
1107     }
1108   else
1109    {
1110      printf("++++++++ real data analysis \n");
1111    }
1112
1113    printf("+++++++ QA hist %d",fqahist);
1114
1115   //---- Geometry
1116   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
1117
1118   //--------Initialize PID
1119   //fPID->SetHasMCData(kFALSE);
1120   fPID->SetHasMCData(fmcData);
1121   if(!fPID->GetNumberOfPIDdetectors())
1122     {
1123       fPID->AddDetector("TPC", 0);
1124       //fPID->AddDetector("EMCAL", 1); <--- apply PID selection in task
1125     }
1126
1127   Double_t params[4];
1128   const char *cutmodel;
1129   cutmodel = "pol0";
1130   params[0] = -1.0; //sigma min
1131   double maxnSig = 3.0;
1132   if(fmcData)
1133     {
1134      params[0] = -5.0; //sigma min
1135      maxnSig = 5.0;
1136     }
1137
1138   for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
1139
1140   fPID->SortDetectors();
1141   fPIDqa = new AliHFEpidQAmanager();
1142   fPIDqa->Initialize(fPID);
1143
1144   //------- fcut --------------
1145   fCuts = new AliHFEcuts();
1146   fCuts->CreateStandardCuts();
1147   //fCuts->SetMinNClustersTPC(100);
1148   fCuts->SetMinNClustersTPC(90);
1149   fCuts->SetMinRatioTPCclusters(0.6);
1150   fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
1151   //fCuts->SetMinNClustersITS(3);
1152   fCuts->SetMinNClustersITS(2);
1153   fCuts->SetProductionVertex(0,50,0,50);
1154   fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
1155   fCuts->SetCheckITSLayerStatus(kFALSE);
1156   fCuts->SetVertexRange(10.);
1157   fCuts->SetTOFPIDStep(kFALSE);
1158   //fCuts->SetPtRange(2, 50);
1159   fCuts->SetPtRange(0.1, 50);
1160   //fCuts->SetMaxImpactParam(3.,3.);
1161   fCuts->SetMaxImpactParam(2.4,3.2); // standard in 2011
1162
1163   //--------Initialize correction Framework and Cuts
1164   fCFM = new AliCFManager;
1165   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
1166   fCFM->SetNStepParticle(kNcutSteps);
1167   for(Int_t istep = 0; istep < kNcutSteps; istep++)
1168     fCFM->SetParticleCutsList(istep, NULL);
1169
1170   if(!fCuts){
1171     AliWarning("Cuts not available. Default cuts will be used");
1172     fCuts = new AliHFEcuts;
1173     fCuts->CreateStandardCuts();
1174   }
1175   fCuts->Initialize(fCFM);
1176
1177   //---------Output Tlist
1178   fOutputList = new TList();
1179   fOutputList->SetOwner();
1180   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
1181
1182   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
1183   fOutputList->Add(fNoEvents);
1184
1185   Int_t binsE[5] =    {250, 100,  1000, 200,   10};
1186   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5};
1187   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5};
1188   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
1189   //if(fqahist==1)fOutputList->Add(fEMCAccE);
1190
1191   hEMCAccE = new TH2F("hEMCAccE","Cluster Energy",200,0,100,100,0,20);
1192   fOutputList->Add(hEMCAccE);
1193
1194   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
1195   fOutputList->Add(fTrkpt);
1196
1197   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
1198   fOutputList->Add(fTrackPtBefTrkCuts);
1199
1200   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
1201   fOutputList->Add(fTrackPtAftTrkCuts);
1202
1203   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
1204   fOutputList->Add(fTPCnsigma);
1205
1206   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
1207   fOutputList->Add(fTrkEovPBef);
1208
1209   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
1210   fOutputList->Add(fTrkEovPAft);
1211
1212   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
1213   fOutputList->Add(fdEdxBef);
1214
1215   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
1216   fOutputList->Add(fdEdxAft);
1217
1218   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
1219   fOutputList->Add(fIncpT);
1220
1221   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1222   fOutputList->Add(fIncpTM20);
1223
1224   Int_t nBinspho[9] =  { 10,  30,   600, 60,   50,    4,  40,   8,  30};
1225   Double_t minpho[9] = {  0.,  0., -0.1, 40,  0, -0.5,   0,-1.5,   0};
1226   Double_t maxpho[9] = {100., 30.,  0.5, 100,   1,  3.5,   2, 6.5,  30};
1227
1228   fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho);
1229   if(fqahist==1)fOutputList->Add(fInvmassLS);
1230
1231   fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho);
1232   if(fqahist==1)fOutputList->Add(fInvmassULS);
1233
1234   fInvmassLSmc = new THnSparseD("fInvmassLSmc", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho);
1235   if(fqahist==1)fOutputList->Add(fInvmassLSmc);
1236
1237   fInvmassULSmc = new THnSparseD("fInvmassULSmc", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho);
1238   if(fqahist==1)fOutputList->Add(fInvmassULSmc);
1239
1240   fInvmassLSreco = new TH2D("fInvmassLSreco", "Inv mass of LS (e,e) reco; cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1241   fInvmassLSreco->Sumw2();
1242   fOutputList->Add(fInvmassLSreco);
1243
1244   fInvmassULSreco = new TH2D("fInvmassULSreco", "Inv mass of ULS (e,e) reco; cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1245   fInvmassULSreco->Sumw2();
1246   fOutputList->Add(fInvmassULSreco);
1247
1248   fInvmassLSmc0 = new TH2D("fInvmassLSmc0", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1249   fInvmassLSmc0->Sumw2();
1250   fOutputList->Add(fInvmassLSmc0);
1251
1252   fInvmassLSmc1 = new TH2D("fInvmassLSmc1", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1253   fInvmassLSmc1->Sumw2();
1254   fOutputList->Add(fInvmassLSmc1);
1255
1256   fInvmassLSmc2 = new TH2D("fInvmassLSmc2", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1257   fInvmassLSmc2->Sumw2();
1258   fOutputList->Add(fInvmassLSmc2);
1259
1260   fInvmassLSmc3 = new TH2D("fInvmassLSmc3", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1261   fInvmassLSmc3->Sumw2();
1262   fOutputList->Add(fInvmassLSmc3);
1263
1264   fInvmassULSmc0 = new TH2D("fInvmassULSmc0", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1265   fInvmassULSmc0->Sumw2();
1266   fOutputList->Add(fInvmassULSmc0);
1267
1268   fInvmassULSmc1 = new TH2D("fInvmassULSmc1", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1269   fInvmassULSmc1->Sumw2();
1270   fOutputList->Add(fInvmassULSmc1);
1271
1272   fInvmassULSmc2 = new TH2D("fInvmassULSmc2", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1273   fInvmassULSmc2->Sumw2();
1274   fOutputList->Add(fInvmassULSmc2);
1275
1276   fInvmassULSmc3 = new TH2D("fInvmassULSmc3", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 );
1277   fInvmassULSmc3->Sumw2();
1278   fOutputList->Add(fInvmassULSmc3);
1279
1280   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1281   fOutputList->Add(fOpeningAngleLS);
1282
1283   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1284   fOutputList->Add(fOpeningAngleULS);
1285
1286   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
1287   fOutputList->Add(fPhotoElecPt);
1288
1289   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
1290   fOutputList->Add(fPhoElecPt);
1291
1292   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
1293   fOutputList->Add(fPhoElecPtM20);
1294
1295   fPhoElecPtM20Mass = new TH2F("fPhoElecPtM20Mass", "Pho-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1296   fPhoElecPtM20Mass->Sumw2();
1297   fOutputList->Add(fPhoElecPtM20Mass);
1298
1299   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
1300   fOutputList->Add(fSameElecPt);
1301
1302   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
1303   fOutputList->Add(fSameElecPtM20);
1304
1305   fSameElecPtM20Mass = new TH2F("fSameElecPtM20Mass", "Same-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1306   fSameElecPtM20Mass->Sumw2();
1307   fOutputList->Add(fSameElecPtM20Mass);
1308
1309   fSemiElecPtM20 = new TH2F("fSemiElecPtM20", "Semi-inclusive electron pt with M20",200,0,100,100,0,50);
1310   fOutputList->Add(fSemiElecPtM20);
1311
1312   fCent = new TH1F("fCent","Centrality",200,0,100) ;
1313   fOutputList->Add(fCent);
1314
1315   // Make common binning
1316   const Double_t kMinP = 0.;
1317   const Double_t kMaxP = 20.;
1318
1319   //+++ 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
1320   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, dca_xy, dca_z, M20, Disp, Centrality, select)
1321   Int_t nBins[16] =  {  100,     7,  60,    20,    90,  100,   25,     60,    60,  100,    40,  10,  250,  100, 100,    8};
1322   Double_t min[16] = {kMinP,  -0.5, 1.0,  -1.0,  -5.0,    0,    0,   -3.0,  -3.0,  0.0,   0.0,   0,    0,    0,  80, -1.5};
1323   Double_t max[16] = {kMaxP,   6.5, 4.0,   1.0,   4.0,  2.0, 0.05,    3.0,   3.0,  1.0,  20.0, 100,  100,  2.0, 180,  6.5};
1324   //fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;clsF;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max);
1325   fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;DCA_xy;DCA_z;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max);
1326   //if(fqahist==1)fOutputList->Add(fEleInfo);
1327
1328   // Make common binning
1329   Int_t nBinsEop[3] =  { 10, 50, 100};
1330   Double_t minEop[3] = {  0,  0,  -5};
1331   Double_t maxEop[3] = {100, 50,   5};
1332   fElenSigma= new THnSparseD("fElenSigma", "Electron nSigma; cent; pT [GeV/c]; nSigma", 3, nBinsEop, minEop, maxEop);
1333   fOutputList->Add(fElenSigma);
1334
1335
1336   //<---  Trigger info
1337   /*
1338   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
1339   fOutputList->Add(fClsEBftTrigCut);
1340
1341   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
1342   fOutputList->Add(fClsEAftTrigCut);
1343
1344   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
1345   fOutputList->Add(fClsEAftTrigCut1);
1346
1347   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
1348   fOutputList->Add(fClsEAftTrigCut2);
1349
1350   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
1351   fOutputList->Add(fClsEAftTrigCut3);
1352
1353   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
1354   fOutputList->Add(fClsEAftTrigCut4);
1355
1356   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1357   fOutputList->Add(fClsETime);
1358
1359   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1360   fOutputList->Add(fClsETime1);
1361
1362   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
1363   fOutputList->Add(fTrigTimes);
1364
1365   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
1366   fOutputList->Add(fCellCheck);
1367   */
1368   //<---------- MC
1369
1370   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1371   fOutputList->Add(fInputHFEMC);
1372
1373   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
1374   fOutputList->Add(fInputAlle);
1375
1376   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1377   fOutputList->Add(fIncpTMChfe);
1378
1379   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
1380   fOutputList->Add(fIncpTMChfeAll);
1381
1382   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1383   fOutputList->Add(fIncpTMCM20hfe);
1384
1385   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
1386   fOutputList->Add(fIncpTMCM20hfeAll);
1387
1388   fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1389   fOutputList->Add(fIncpTMCM20hfeCheck);
1390
1391   Int_t nBinspho2[5] =  { 200, 100,    7,    3, 700};
1392   Double_t minpho2[5] = {  0.,  0., -2.5, -0.5, 0.};
1393   Double_t maxpho2[5] = {100., 50.,  4.5,  2.5, 70.};
1394
1395   fInputHFEMC_weight = new THnSparseD("fInputHFEMC_weight", "MC HFE electron pt",5,nBinspho2,minpho2,maxpho2);
1396   fOutputList->Add(fInputHFEMC_weight);
1397
1398   fIncpTMCM20hfeCheck_weight = new THnSparseD("fIncpTMCM20hfeCheck_weight", "HFE electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1399   fOutputList->Add(fIncpTMCM20hfeCheck_weight);
1400
1401   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
1402   fOutputList->Add(fIncpTMCpho);
1403
1404   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1405   fOutputList->Add(fIncpTMCM20pho);
1406
1407   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1408   fOutputList->Add(fPhoElecPtMC);
1409
1410   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1411   fOutputList->Add(fPhoElecPtMCM20);
1412
1413   fPhoElecPtMCM20Mass = new TH2D("fPhoElecPtMCM20Mass", "MC Pho-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1414   fOutputList->Add(fPhoElecPtMCM20Mass);
1415
1416   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1417   fOutputList->Add(fSameElecPtMC);
1418
1419   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1420   fOutputList->Add(fSameElecPtMCM20);
1421
1422   fSameElecPtMCM20Mass = new TH2D("fSameElecPtMCM20Mass", "MC Same-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1423   fOutputList->Add(fSameElecPtMCM20Mass);
1424
1425   fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1426   fIncpTMCM20pho_pi0e->Sumw2();
1427   fOutputList->Add(fIncpTMCM20pho_pi0e);
1428
1429   fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1430   fPhoElecPtMCM20_pi0e->Sumw2();
1431   fOutputList->Add(fPhoElecPtMCM20_pi0e);
1432
1433   fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1434   fSameElecPtMCM20_pi0e->Sumw2();
1435   fOutputList->Add(fSameElecPtMCM20_pi0e);
1436  //
1437   fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1438   fIncpTMCM20pho_eta->Sumw2();
1439   fOutputList->Add(fIncpTMCM20pho_eta);
1440
1441   fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1442   fPhoElecPtMCM20_eta->Sumw2();
1443   fOutputList->Add(fPhoElecPtMCM20_eta);
1444
1445   fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1446   fSameElecPtMCM20_eta->Sumw2();
1447   fOutputList->Add(fSameElecPtMCM20_eta);
1448   // ------------
1449   fIncpTMCpho_pi0e_TPC = new THnSparseD("fIncpTMCpho_pi0e_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1450   fIncpTMCpho_pi0e_TPC->Sumw2();
1451   fOutputList->Add(fIncpTMCpho_pi0e_TPC);
1452
1453   fPhoElecPtMC_pi0e_TPC = new THnSparseD("fPhoElecPtMC_pi0e_TPC", "MC Pho-inclusive electron pt with  pi0->e",5,nBinspho2,minpho2,maxpho2);
1454   fPhoElecPtMC_pi0e_TPC->Sumw2();
1455   fOutputList->Add(fPhoElecPtMC_pi0e_TPC);
1456
1457   fSameElecPtMC_pi0e_TPC = new THnSparseD("fSameElecPtMC_pi0e_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1458   fSameElecPtMC_pi0e_TPC->Sumw2();
1459   fOutputList->Add(fSameElecPtMC_pi0e_TPC);
1460
1461   fIncpTMCpho_eta_TPC = new THnSparseD("fIncpTMCpho_eta_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1462   fIncpTMCpho_eta_TPC->Sumw2();
1463   fOutputList->Add(fIncpTMCpho_eta_TPC);
1464
1465   fPhoElecPtMC_eta_TPC = new THnSparseD("fPhoElecPtMC_eta_TPC", "MC Pho-inclusive electron pt with  pi0->e",5,nBinspho2,minpho2,maxpho2);
1466   fPhoElecPtMC_eta_TPC->Sumw2();
1467   fOutputList->Add(fPhoElecPtMC_eta_TPC);
1468
1469   fSameElecPtMC_eta_TPC = new THnSparseD("fSameElecPtMC_eta_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1470   fSameElecPtMC_eta_TPC->Sumw2();
1471   fOutputList->Add(fSameElecPtMC_eta_TPC);
1472
1473
1474   //-------------
1475
1476   CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1477   fOutputList->Add(CheckNclust);
1478
1479   CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1480   fOutputList->Add(CheckNits);
1481
1482   CheckDCA = new TH2D("CheckDCA","DCA check",200,-5,5,200,-5,5);
1483   fOutputList->Add(CheckDCA);
1484   /*
1485   Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1486   fOutputList->Add(Hpi0pTcheck);
1487
1488   HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1489   fOutputList->Add(HETApTcheck);
1490   */
1491
1492   Int_t nBinspho3[3] =  { 200, 100, 3};
1493   Double_t minpho3[3] = {  0.,  0., -0.5};
1494   Double_t maxpho3[3] = {100., 50., 2.5};
1495
1496   Hpi0pTcheck = new THnSparseD("Hpi0pTcheck","Pi0 pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1497   fOutputList->Add(Hpi0pTcheck);
1498
1499   HETApTcheck = new THnSparseD("HETApTcheck","Eta pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1500   fOutputList->Add(HETApTcheck);
1501   //--
1502   HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1503   fOutputList->Add(HphopTcheck);
1504   //
1505   HDpTcheck = new TH2D("HDpTcheck","D pT from Hijing",100,0,50,3,-0.5,2.5);
1506   fOutputList->Add(HDpTcheck);
1507
1508   HBpTcheck = new TH2D("HBpTcheck","B pT from Hijing",100,0,50,3,-0.5,2.5);
1509   fOutputList->Add(HBpTcheck);
1510   //
1511
1512   fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1513   fOutputList->Add(fpTCheck);
1514
1515   fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1516   fOutputList->Add(fMomDtoE);
1517
1518   fLabelCheck = new TH2D("fLabelCheck","MC label",50,0,50,5,-1.5,3.5);
1519   fOutputList->Add(fLabelCheck);
1520
1521   fgeoFake = new TH2D("fgeoFake","Label==0 eta and phi",628,0,6.28,200,-1,1);
1522   fOutputList->Add(fgeoFake);
1523
1524   fFakeTrk0 = new TH2D("fFakeTrk0","fake trakcs",10,0,100,20,0,20);
1525   fOutputList->Add(fFakeTrk0);
1526
1527   fFakeTrk1 = new TH2D("fFakeTrk1","true all e a.f. eID",10,0,100,20,0,20);
1528   fOutputList->Add(fFakeTrk1);
1529
1530   ftimingEle = new TH2D("ftimingEle","electron TOF",100,0,20,100,1e-7,1e-6);
1531   fOutputList->Add(ftimingEle);
1532
1533   // eta correction
1534   // note: parameters 01/31new.TPCnSigmaEtaDep
1535   // 70-90 delta_eta = 0.2
1536
1537   double etaval[12] = {-0.55,-0.45,-0.35,-0.25,-0.15,-0.05,0.05,0.15,0.25,0.35,0.45,0.55};
1538   double corr0[12]= {-0.569177,-0.528844,-0.391979,-0.165494,0.0283495,0.156171,0.266353,0.13103,-0.0250842,-0.274089,-0.45481,-0.536291}; // 0-10 (done)
1539   double corr1[12]= {-0.404742,-0.278953,-0.218069,0.00139927,0.191412,0.354403,0.524594,0.341778,0.244199,-0.112146,-0.160692,-0.352832}; // 10-20 (done)
1540   double corr2[12] = {-0.306007,-0.16821,-0.0248635,0.202233,0.447051,0.497197,0.712251,0.433482,0.337907,0.168426,-0.0693229,-0.0728351}; // 20-30 (done)
1541   double corr3[12] = {-0.13884,-0.0503553,0.104403,0.389773,0.50697,0.539048,0.751642,0.655636,0.518563,0.308156,0.0361159,-0.0491439}; // 30-40 (done)
1542   double corr4[12] = {-0.0319431,0.0808711,0.208774,0.443217,0.557762,0.61453,0.889519,0.808282,0.620394,0.267092,0.15241,-0.0458664}; // 40-50 (done)
1543   double corr5[12] = {-0.130625,0.0189124,0.190344,0.467431,0.546353,0.672251,0.731541,0.802101,0.437108,0.294081,0.193682,0.159074}; // 50-70(done)
1544   double corr6[12] = {0.0600197,0.0600197,0.358366,0.358366,0.973734,0.973734,0.759812,0.759812,0.667861,0.667861,0.415635,0.415635}; // 70-90(done)
1545
1546   fnSigEtaCorr[0] = new TGraphErrors(12,etaval,corr0); // 0-10
1547   fnSigEtaCorr[1] = new TGraphErrors(12,etaval,corr1); // 10-20
1548   fnSigEtaCorr[2] = new TGraphErrors(12,etaval,corr2); // 20-30
1549   fnSigEtaCorr[3] = new TGraphErrors(12,etaval,corr3); // 30-40
1550   fnSigEtaCorr[4] = new TGraphErrors(12,etaval,corr4); // 40-50
1551   fnSigEtaCorr[5] = new TGraphErrors(12,etaval,corr5); // 50-70
1552   fnSigEtaCorr[6] = new TGraphErrors(12,etaval,corr6); // 70-90
1553
1554   fIncMaxE = new TH2D("fIncMaxE","Inc",10,0,100,10,0,100);
1555   fOutputList->Add(fIncMaxE);
1556
1557   fIncReco = new TH2D("fIncReco","Inc",10,0,100,100,0,500);
1558   fOutputList->Add(fIncReco);
1559
1560   fPhoReco = new TH2D("fPhoReco","Pho",10,0,100,100,0,500);
1561   fOutputList->Add(fPhoReco);
1562
1563   fSamReco = new TH2D("fSamReco","Same",10,0,100,100,0,500);
1564   fOutputList->Add(fSamReco);
1565
1566   fIncRecoMaxE = new TH2D("fIncRecoMaxE","Inc",10,0,100,100,0,500);
1567   fOutputList->Add(fIncRecoMaxE);
1568
1569   fPhoRecoMaxE = new TH2D("fPhoRecoMaxE","Pho",10,0,100,100,0,500);
1570   fOutputList->Add(fPhoRecoMaxE);
1571
1572   fSamRecoMaxE = new TH2D("fSamRecoMaxE","Same",10,0,100,100,0,500);
1573   fOutputList->Add(fSamRecoMaxE);
1574
1575   fPhoVertexReco_TPC = new TH2D("fPhoVertexReco_TPC","photon production Vertex mass selection TPC",40,0,20,250,0,50);
1576   fPhoVertexReco_TPC->Sumw2();
1577   fOutputList->Add(fPhoVertexReco_TPC);
1578
1579   fPhoVertexReco_TPC_Invmass = new TH2D("fPhoVertexReco_TPC_Invmass","photon production Vertex mass selection TPC Invmass",40,0,20,250,0,50);
1580   fPhoVertexReco_TPC_Invmass->Sumw2();
1581   fOutputList->Add(fPhoVertexReco_TPC_Invmass);
1582
1583   fPhoVertexReco_EMCal = new TH2D("fPhoVertexReco_EMCal","photon production Vertex mass selection",40,0,20,250,0,50);
1584   fPhoVertexReco_EMCal->Sumw2();
1585   fOutputList->Add(fPhoVertexReco_EMCal);
1586
1587   fPhoVertexReco_Invmass = new TH2D("fPhoVertexReco_Invmass","photon production Vertex mass selection",40,0,20,250,0,50);
1588   fPhoVertexReco_Invmass->Sumw2();
1589   fOutputList->Add(fPhoVertexReco_Invmass);
1590
1591   fPhoVertexReco_step0= new TH2D("fPhoVertexReco_step0","photon production Vertex mass selection",40,0,20,250,0,50);
1592   fPhoVertexReco_step0->Sumw2();
1593   fOutputList->Add(fPhoVertexReco_step0);
1594
1595   fPhoVertexReco_step1= new TH2D("fPhoVertexReco_step1","photon production Vertex mass selection",40,0,20,250,0,50);
1596   fPhoVertexReco_step1->Sumw2();
1597   fOutputList->Add(fPhoVertexReco_step1);
1598
1599   fMatchV0_0 = new TH1D("fMatchV0_0","V0 match",100,0,20);
1600   fOutputList->Add(fMatchV0_0);
1601
1602   fMatchV0_1 = new TH1D("fMatchV0_1","V0 match",100,0,20);
1603   fOutputList->Add(fMatchV0_1);
1604
1605   fMatchMC_0 = new TH1D("fMatchMC_0","MC match",100,0,20);
1606   fOutputList->Add(fMatchMC_0);
1607
1608   fMatchMC_1 = new TH1D("fMatchMC_1","MC match",100,0,20);
1609   fOutputList->Add(fMatchMC_1);
1610
1611   fpair = new TH2D("fpair","pair of associate",100,0,20,21,-10.5,10.5);
1612   fOutputList->Add(fpair);
1613
1614   fFakeRejection0 = new TH2D("fFakeRejection0","TPC PID",2,-0.5,1.5,100,0,20);
1615   fFakeRejection0->Sumw2();
1616   fOutputList->Add(fFakeRejection0);
1617
1618   fFakeRejection1 = new TH2D("fFakeRejection1","TPC PID + Tr match",2,-0.5,1.5,100,0,20);
1619   fFakeRejection1->Sumw2();
1620   fOutputList->Add(fFakeRejection1);
1621
1622   fFakeRejection2 = new TH2D("fFakeRejection2","TPC PID + Tr match + E/p",2,-0.5,1.5,100,0,20);
1623   fFakeRejection2->Sumw2();
1624   fOutputList->Add(fFakeRejection2);
1625
1626   EopFake = new TH2D("EopFake","negative track Eop",200,0,20,200,0,2);
1627   EopFake->Sumw2();
1628   fOutputList->Add(EopFake);
1629
1630   EopTrue = new TH2D("EopTrue","true track Eop",200,0,20,200,0,2);
1631   EopTrue->Sumw2();
1632   fOutputList->Add(EopTrue);
1633
1634   MatchFake = new TH2D("MatchFake","negative track Match",200,0,20,100,0,0.05);
1635   MatchFake->Sumw2();
1636   fOutputList->Add(MatchFake);
1637
1638   MatchTrue = new TH2D("MatchTrue","true track Match",200,0,20,100,0,05);
1639   MatchTrue->Sumw2();
1640   fOutputList->Add(MatchTrue);
1641
1642   PostData(1,fOutputList);
1643 }
1644
1645 //________________________________________________________________________
1646 void AliAnalysisTaskHFECal::Terminate(Option_t *)
1647 {
1648   // Info("Terminate");
1649         AliAnalysisTaskSE::Terminate();
1650 }
1651
1652 //________________________________________________________________________
1653 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1654 {
1655   // Check single track cuts for a given cut step
1656   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1657   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1658   return kTRUE;
1659 }
1660 //_________________________________________
1661 void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonic, Bool_t &fFlagConvinat, Double_t nSig, Double_t shower, Double_t ep, Double_t mce, Double_t w, Int_t ibgevent, Bool_t tagpi0, Bool_t tageta, Int_t iCal)
1662 {
1663   //Identify non-heavy flavour electrons using Invariant mass method
1664
1665   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1666   fTrackCuts->SetRequireTPCRefit(kTRUE);
1667   fTrackCuts->SetRequireITSRefit(kTRUE);
1668   fTrackCuts->SetEtaRange(-0.9,0.9);
1669   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1670   fTrackCuts->SetMinNClustersTPC(90);
1671
1672   //const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1673   Double_t bfield = fESD->GetMagneticField();
1674   Double_t emass = 0.51*0.001; // (0.51 MeV)
1675
1676   double ptEle = track->Pt();  //add
1677   if(ibgevent==0 && w > 0.0)
1678      {
1679       fpTCheck->Fill(ptEle,w);
1680      }
1681
1682   Bool_t flagPhotonicElec = kFALSE;
1683   Bool_t flagConvinatElec = kFALSE;
1684
1685   int p1 = 0;
1686   if(mce==3)
1687      {
1688        Int_t label = TMath::Abs(track->GetLabel());
1689        TParticle* particle = stack->Particle(label);
1690        p1 = particle->GetFirstMother();
1691      }
1692
1693   int numULS = 0;
1694   int numLS = 0;
1695
1696   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1697     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1698     if (!trackAsso) {
1699       printf("ERROR: Could not receive track %d\n", jTracks);
1700       continue;
1701     }
1702     if(itrack==jTracks)continue;
1703     int jbgevent = 0;
1704
1705     int p2 = 0;
1706     if(mce==3)
1707     {
1708       Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1709       TParticle* particle2 = stack->Particle(label2);
1710       Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1711       if(MChijing_ass)jbgevent =1;
1712       if(particle2->GetFirstMother()>-1)
1713          p2 = particle2->GetFirstMother();
1714     }
1715
1716     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1717     Double_t mass=999., width = -999;
1718     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1719
1720     //ptPrim = track->Pt();
1721     ptPrim = ptEle;
1722
1723     dEdxAsso = trackAsso->GetTPCsignal();
1724     ptAsso = trackAsso->Pt();
1725     Int_t chargeAsso = trackAsso->Charge();
1726     Int_t charge = track->Charge();
1727
1728
1729     if(ptAsso <fMimpTassCut) continue;
1730     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1731     double fTPCnSigmaAss = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
1732     //cout << "fTPCnSigmaAss = " << fTPCnSigmaAss << endl;
1733     //cout << "fTPCnSigmaAss Cut = " << fMimNsigassCut << endl;
1734     if(fTPCnSigmaAss <fMimNsigassCut || fTPCnSigmaAss>5) continue;
1735     //cout << "fTPCnSigmaAss a.f. cut = " << fTPCnSigmaAss << endl;
1736
1737     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1738     if(charge>0) fPDGe1 = -11;
1739     if(chargeAsso>0) fPDGe2 = -11;
1740
1741     //printf("chargeAsso = %d\n",chargeAsso);
1742     //printf("charge = %d\n",charge);
1743     if(charge == chargeAsso) fFlagLS = kTRUE;
1744     if(charge != chargeAsso) fFlagULS = kTRUE;
1745
1746     //printf("fFlagLS = %d\n",fFlagLS);
1747     //printf("fFlagULS = %d\n",fFlagULS);
1748     //printf("\n");
1749
1750     AliKFParticle::SetField(bfield);
1751     AliKFParticle ge1(*track, fPDGe1);
1752     AliKFParticle ge2(*trackAsso, fPDGe2);
1753
1754     if(fSetMassNonlinear)
1755       {
1756        ge1.SetNonlinearMassConstraint(emass);
1757        ge2.SetNonlinearMassConstraint(emass);
1758       }
1759
1760     AliKFParticle recg(ge1, ge2);
1761
1762     /*
1763     // vertex
1764     AliKFVertex primV(*pVtx);
1765     primV += recg;
1766     primV -= ge1;
1767     primV -= ge2;
1768     recg.SetProductionVertex(primV);
1769     */
1770
1771     // check chi2
1772     if(recg.GetNDF()<1) continue;
1773     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1774     Double_t chi2cut = 3.0;
1775
1776     // mass.
1777     if(fSetMassConstraint)
1778       {
1779        recg.SetMassConstraint(0,0.0001);
1780        chi2cut = 30.0;
1781       }
1782     recg.GetMass(mass,width);
1783
1784     if(fSetMassWidthCut && width>1e10)continue;
1785
1786         // angle
1787     openingAngle = ge1.GetAngle(ge2);
1788     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1789     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1790
1791     double ishower = 0;
1792     if(shower>0.0 && shower<0.3)ishower = 1;
1793
1794     if(iCal==1)
1795       {
1796         double phoinfo[9];
1797         phoinfo[0] = cent;
1798         phoinfo[1] = ptPrim;
1799         phoinfo[2] = mass;
1800         phoinfo[3] = nSig;
1801         phoinfo[4] = openingAngle;
1802         phoinfo[5] = ishower;
1803         phoinfo[6] = ep;
1804         phoinfo[7] = mce;
1805         phoinfo[8] = ptAsso;
1806
1807         if(fFlagLS) fInvmassLS->Fill(phoinfo);
1808         if(fFlagULS) fInvmassULS->Fill(phoinfo);
1809         if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1810         if(fFlagULS && ibgevent==0 && jbgevent==0)
1811         {
1812               fInvmassULSmc->Fill(phoinfo,w);
1813         }
1814        //printf("fInvmassCut %f\n",fInvmassCut);
1815        //printf("openingAngle %f\n",fOpeningAngleCut);
1816       }
1817
1818     // angle cut
1819     if(openingAngle > fOpeningAngleCut) continue;
1820     // chi2 cut
1821     //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue;
1822     if(chi2recg>chi2cut) continue;
1823
1824     if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass);
1825     if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass);
1826
1827     // check double count
1828     if(mass<fInvmassCut && fFlagULS)numULS++;
1829     if(mass<fInvmassCut && fFlagLS)numLS++;
1830
1831     // for real data
1832     //printf("mce =%f\n",mce);
1833
1834     if(mce<-0.5) // mce==-1. is real
1835     {
1836             //printf("Real data\n");
1837
1838             // count
1839             if(mass<fInvmassCut && fFlagULS && shower>0.0 && shower<0.3 && iCal==1)fPhoElecPtM20Mass->Fill(cent,ptPrim);
1840             if(mass<fInvmassCut && fFlagLS && shower>0.0 && shower<0.3 && iCal==1)fSameElecPtM20Mass->Fill(cent,ptPrim);
1841
1842             // flag
1843             if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1844                     flagPhotonicElec = kTRUE;
1845             }
1846             if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1847                     flagConvinatElec = kTRUE;
1848             }
1849     }
1850     // for MC data
1851     else
1852     {
1853             // count
1854             if(mass<fInvmassCut && fFlagULS && shower>0.0 && shower<0.3 && mce>2.1 && iCal==1)fPhoElecPtMCM20Mass->Fill(cent,ptPrim,w);
1855             if(mass<fInvmassCut && fFlagLS && shower>0.0 && shower<0.3 && mce>2.1 && iCal==1)fSameElecPtMCM20Mass->Fill(cent,ptPrim,w);
1856
1857             //printf("MC data\n");
1858
1859             if(w>0.0)
1860             {
1861                     //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1862                     if(iCal==1)
1863                       {
1864                        if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass);
1865                        if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass);
1866                        if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass);
1867                        if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass);
1868                        if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass);
1869                        if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass);
1870                        if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass);
1871                        if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass);
1872                       }
1873             }
1874
1875             if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1876                     //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ //<--- only MC train (55,56) v5-03-68-AN , 69 & v5-05-70-AN (till 74)
1877                     flagPhotonicElec = kTRUE;
1878             }
1879             if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1880                     flagConvinatElec = kTRUE;
1881             }
1882     }
1883
1884   } // end of associate loop
1885
1886   if(numULS>0 || numLS>0)
1887     {
1888      int numPair = numULS-numLS;
1889      if(iCal==1)fpair->Fill(ptEle,numPair);
1890     }
1891
1892   fFlagPhotonic = flagPhotonicElec;
1893   fFlagConvinat = flagConvinatElec;
1894
1895 }
1896
1897 //-------------------------------------------
1898
1899 void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1900 {
1901  //int label = 99999;
1902  //int pid = 99999;
1903
1904  if(part->GetFirstMother()>-1)
1905    {
1906     label = part->GetFirstMother();
1907     pid = stack->Particle(label)->GetPdgCode();
1908    }
1909    //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
1910 }
1911
1912 double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1913 {
1914         double weight = 1.0;
1915
1916         if(mcPi0pT>0.0 && mcPi0pT<5.0)
1917         {
1918                 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1919         }
1920         else
1921         {
1922                 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1923         }
1924   return weight;
1925 }
1926
1927 double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1928 {
1929   double weight = 1.0;
1930
1931   weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1932   return weight;
1933 }
1934
1935
1936 //_________________________________________
1937 void AliAnalysisTaskHFECal::FindTriggerClusters()
1938 {
1939   //cout << "finding trigger patch" << endl;
1940   // constants
1941   const int nModuleCols = 2;
1942   const int nModuleRows = 5;
1943   const int nColsFeeModule = 48;
1944   const int nRowsFeeModule = 24;
1945   const int nColsFaltroModule = 24;
1946   const int nRowsFaltroModule = 12;
1947   //const int faltroWidthMax = 20;
1948
1949   // part 1, trigger extraction -------------------------------------
1950   Int_t globCol, globRow;
1951   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1952   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1953
1954   //Int_t trigtimes[faltroWidthMax];
1955   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1956   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1957   //Double_t fTrigCutLow = 6;
1958   //Double_t fTrigCutHigh = 10;
1959   Double_t fTimeCutLow =  469e-09;
1960   Double_t fTimeCutHigh = 715e-09;
1961
1962   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1963
1964   // erase trigger maps
1965   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1966   {
1967     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
1968     {
1969       ftriggersCut[i][j] = 0;
1970       ftriggers[i][j] = 0;
1971       ftriggersTime[i][j] = 0;
1972     }
1973   }
1974
1975   Int_t iglobCol=0, iglobRow=0;
1976   // go through triggers
1977   if( fCaloTrigger->GetEntries() > 0 )
1978   {
1979     // needs reset
1980     fCaloTrigger->Reset();
1981     while( fCaloTrigger->Next() )
1982     {
1983       fCaloTrigger->GetPosition( globCol, globRow );
1984       fCaloTrigger->GetNL0Times( ntimes );
1985       /*
1986       // no L0s
1987       if( ntimes < 1 )   continue;
1988       // get precise timings
1989       fCaloTrigger->GetL0Times( trigtimes );
1990       trigInCut = 0;
1991       for(Int_t i = 0; i < ntimes; i++ )
1992       {
1993         // save the first trigger time in channel
1994         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
1995           triggersTime[globCol][globRow] = trigtimes[i];
1996         //printf("trigger times: %d\n",trigtimes[i]);
1997         // check if in cut
1998         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
1999           trigInCut = 1;
2000
2001         fTrigTimes->Fill(trigtimes[i]);
2002       }
2003       */
2004
2005       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
2006       Int_t bit = 0;
2007       fCaloTrigger->GetTriggerBits(bit);
2008       //cout << "bit = " << bit << endl;
2009
2010       Int_t ts = 0;
2011       fCaloTrigger->GetL1TimeSum(ts);
2012       //cout << "ts = " << ts << endl;
2013       if (ts > 0)ftriggers[globCol][globRow] = 1;
2014       // number of triggered channels in event
2015       nTrigChannel++;
2016       // ... inside cut
2017       if(ts>0 && (bit >> 6 & 0x1))
2018       {
2019         iglobCol = globCol;
2020         iglobRow = globRow;
2021         nTrigChannelCut++;
2022         //cout << "ts cut = " << ts << endl;
2023         //cout << "globCol = " << globCol << endl;
2024         //cout << "globRow = " << globRow << endl;
2025         ftriggersCut[globCol][globRow] = 1;
2026       }
2027
2028     } // calo trigger entries
2029   } // has calo trigger entries
2030
2031   // part 2 go through the clusters here -----------------------------------
2032   //cout << " part 2 go through the clusters here ----------------------------------- " << endl;
2033   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
2034   Short_t cellAddr, nSACell;
2035   Int_t mclabel;
2036   //Int_t nSACell, iSACell, mclabel;
2037   Int_t iSACell;
2038   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
2039   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
2040
2041   TRefArray *fCaloClusters = new TRefArray();
2042   fESD->GetEMCALClusters( fCaloClusters );
2043   nCluster = fCaloClusters->GetEntries();
2044
2045
2046   // save all cells times if there are clusters
2047   if( nCluster > 0 ){
2048     // erase time map
2049     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
2050       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
2051         cellTime[i][j] = 0.;
2052         cellEnergy[i][j] = 0.;
2053       }
2054     }
2055
2056     // get cells
2057     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
2058     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
2059     nSACell = fCaloCells->GetNumberOfCells();
2060     for(iSACell = 0; iSACell < nSACell; iSACell++ ){
2061       // get the cell info *fCal
2062       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
2063
2064       // get cell position
2065       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
2066       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2067
2068       // convert co global phi eta
2069       gphi = iphi + nRowsFeeModule*(nSupMod/2);
2070       geta = ieta + nColsFeeModule*(nSupMod%2);
2071
2072       // save cell time and energy
2073       cellTime[geta][gphi] = cellTimeT;
2074       cellEnergy[geta][gphi] = cellAmp;
2075
2076     }
2077   }
2078
2079   Int_t nClusterTrig, nClusterTrigCut;
2080   UShort_t *cellAddrs;
2081   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
2082   Float_t clsPos[3] = {0.,0.,0.};
2083
2084   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
2085   {
2086     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2087     if(!cluster || !cluster->IsEMCAL()) continue;
2088
2089     // get cluster cells
2090     nCell = cluster->GetNCells();
2091
2092     // get cluster energy
2093     clsE = cluster->E();
2094
2095     // get cluster position
2096     cluster->GetPosition(clsPos);
2097     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2098     clsEta = clsPosVec.Eta();
2099     clsPhi = clsPosVec.Phi();
2100
2101     // get the cell addresses
2102     cellAddrs = cluster->GetCellsAbsId();
2103
2104     // check if the cluster contains cell, that was marked as triggered
2105     nClusterTrig = 0;
2106     nClusterTrigCut = 0;
2107
2108     // loop the cells to check, if cluser in acceptance
2109     // any cluster with a cell outside acceptance is not considered
2110     for( iCell = 0; iCell < nCell; iCell++ )
2111     {
2112      // check hot cell
2113      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
2114
2115       // get cell position
2116       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
2117       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2118
2119       // convert co global phi eta
2120       gphi = iphi + nRowsFeeModule*(nSupMod/2);
2121       geta = ieta + nColsFeeModule*(nSupMod%2);
2122
2123       if( cellTime[geta][gphi] > 0. ){
2124         clusterTime += cellTime[geta][gphi];
2125         gCell++;
2126       }
2127
2128       // get corresponding FALTRO
2129       fphi = gphi / 2;
2130       feta = geta / 2;
2131
2132         //cout << "fphi = " << fphi << endl;
2133         //cout << "feta = " << feta << endl;
2134
2135       // try to match with a triggered
2136       if( ftriggers[feta][fphi]==1)
2137       {  nClusterTrig++;
2138       }
2139       if( ftriggersCut[feta][fphi]==1)
2140       { nClusterTrigCut++;
2141       }
2142
2143       //cout << "nClusterTrigCut : " << nClusterTrigCut << endl;
2144
2145     } // cells
2146
2147
2148     if( gCell > 0 )
2149       clusterTime = clusterTime / (Double_t)gCell;
2150     // fix the reconstriction code time 100ns jumps
2151     if( fESD->GetBunchCrossNumber() % 4 < 2 )
2152       clusterTime -= 0.0000001;
2153
2154     //fClsETime->Fill(clsE,clusterTime);
2155     //fClsEBftTrigCut->Fill(clsE);
2156
2157     if(nClusterTrig>0){
2158       //fClsETime1->Fill(clsE,clusterTime);
2159     }
2160
2161     if(nClusterTrig>0){
2162       cluster->SetChi2(1);
2163       //fClsEAftTrigCut1->Fill(clsE);
2164     }
2165
2166     if(nClusterTrigCut>0){
2167       cluster->SetChi2(2);
2168       //fClsEAftTrigCut2->Fill(clsE);
2169     }
2170
2171     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
2172     {
2173       cluster->SetChi2(3);
2174       //fClsEAftTrigCut3->Fill(clsE);
2175     }
2176
2177     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
2178     {
2179       // cluster->SetChi2(4);
2180       //fClsEAftTrigCut4->Fill(clsE);
2181     }
2182     if(nClusterTrigCut<1)
2183     {
2184       cluster->SetChi2(0);
2185
2186       //fClsEAftTrigCut->Fill(clsE);
2187     }
2188
2189   } // clusters
2190 }
2191
2192 // <-------- only MC correction
2193 double AliAnalysisTaskHFECal::MCEopMeanCorrection(double pTmc, float central)
2194 {
2195   TF1 *fcorr0 = new TF1("fcorr0","[0]*tanh([1]+[2]*x)");
2196   TF1 *fcorr1 = new TF1("fcorr1","[0]*tanh([1]+[2]*x)");
2197
2198  double shift = 0.0;
2199
2200  if(central>0 && central<=10)
2201    {
2202     fcorr0->SetParameters(1.045,1.288,3.18e-01); //
2203     fcorr1->SetParameters(9.91e-01,3.466,2.344);
2204    }
2205  else if(central>10 && central<=20)
2206    {
2207     fcorr0->SetParameters(1.029,8.254e-01,4.07e-01);
2208     fcorr1->SetParameters(0.975,2.276,1.501e-01);
2209    }
2210  else if(central>20 && central<=30)
2211    {
2212     fcorr0->SetParameters(1.01,8.795e-01,3.904e-01);
2213     fcorr1->SetParameters(9.675e-01,1.654,2.583e-01);
2214    }
2215  else if(central>30 && central<=40)
2216    {
2217     fcorr0->SetParameters(1.00,1.466,2.305e-1);
2218     fcorr1->SetParameters(9.637e-01,1.473,2.754e-01);
2219    }
2220  else if(central>40 && central<=50)
2221    {
2222     fcorr0->SetParameters(1.00,1.422,1.518e-01);
2223     fcorr1->SetParameters(9.59e-01,1.421,2.931e-01);
2224    }
2225
2226  else if(central>50 && central<=70)
2227    {
2228     fcorr0->SetParameters(0.989,2.495,2.167);
2229     fcorr1->SetParameters(0.961,1.734,1.438e-01);
2230    }
2231  else if(central>70 && central<=100)
2232    {
2233     fcorr0->SetParameters(0.981,-3.373,3.93327);
2234     fcorr1->SetParameters(9.574e-01,1.698,1.58e-01);
2235    }
2236
2237
2238  shift = fcorr0->Eval(pTmc)-fcorr1->Eval(pTmc);
2239
2240  fcorr0->Delete();
2241  fcorr1->Delete();
2242
2243  return shift;
2244 }
2245
2246 // <-------- only Data correction
2247 double AliAnalysisTaskHFECal::NsigmaCorrection(double tmpeta, float central)
2248 {
2249  int icent = 0;
2250
2251  if(central>=0 && central<10)
2252    {
2253     icent = 0;
2254    }
2255  else if(central>=10 && central<20)
2256   {
2257    icent = 1;
2258   }
2259  else if(central>=20 && central<30)
2260   {
2261    icent = 2;
2262   }
2263  else if(central>=30 && central<40)
2264   {
2265    icent = 3;
2266   }
2267  else if(central>=40 && central<50)
2268   {
2269    icent = 4;
2270   }
2271  else if(central>=50 && central<70)
2272   {
2273    icent = 5;
2274   }
2275  else
2276   {
2277    icent = 6;
2278   }
2279
2280  double shift = fnSigEtaCorr[icent]->Eval(tmpeta);
2281
2282  //cout << "eta correction"<< endl;
2283  //cout << "cent = "<< central<< endl;
2284  //cout << "icent = "<< icent << endl;
2285  //cout << "shift = "<< shift << endl;
2286
2287  return shift;
2288
2289 }
2290
2291