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