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