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); // default
1181   fCuts->SetCutITSpixel(AliHFEextraCuts::kBoth);
1182   fCuts->SetCheckITSLayerStatus(kFALSE);
1183   fCuts->SetVertexRange(10.);
1184   fCuts->SetTOFPIDStep(kFALSE);
1185   //fCuts->SetPtRange(2, 50);
1186   fCuts->SetPtRange(0.1, 50);
1187   //fCuts->SetMaxImpactParam(3.,3.);
1188   fCuts->SetMaxImpactParam(2.4,3.2); // standard in 2011
1189
1190   //--------Initialize correction Framework and Cuts
1191   fCFM = new AliCFManager;
1192   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
1193   fCFM->SetNStepParticle(kNcutSteps);
1194   for(Int_t istep = 0; istep < kNcutSteps; istep++)
1195     fCFM->SetParticleCutsList(istep, NULL);
1196
1197   if(!fCuts){
1198     AliWarning("Cuts not available. Default cuts will be used");
1199     fCuts = new AliHFEcuts;
1200     fCuts->CreateStandardCuts();
1201   }
1202   fCuts->Initialize(fCFM);
1203
1204   //---------Output Tlist
1205   fOutputList = new TList();
1206   fOutputList->SetOwner();
1207   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
1208
1209   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
1210   fOutputList->Add(fNoEvents);
1211
1212   Int_t binsE[5] =    {250, 100,  1000, 200,   10};
1213   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5};
1214   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5};
1215   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
1216   //if(fqahist==1)fOutputList->Add(fEMCAccE);
1217
1218   hEMCAccE = new TH2F("hEMCAccE","Cluster Energy",200,0,100,100,0,20);
1219   fOutputList->Add(hEMCAccE);
1220
1221   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
1222   fOutputList->Add(fTrkpt);
1223
1224   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
1225   fOutputList->Add(fTrackPtBefTrkCuts);
1226
1227   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
1228   fOutputList->Add(fTrackPtAftTrkCuts);
1229
1230   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
1231   fOutputList->Add(fTPCnsigma);
1232
1233   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
1234   fOutputList->Add(fTrkEovPBef);
1235
1236   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
1237   fOutputList->Add(fTrkEovPAft);
1238
1239   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
1240   fOutputList->Add(fdEdxBef);
1241
1242   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
1243   fOutputList->Add(fdEdxAft);
1244
1245   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
1246   fOutputList->Add(fIncpT);
1247
1248   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1249   fOutputList->Add(fIncpTM20);
1250
1251   Int_t nBinspho[9] =  { 10,  30,   600, 60,   50,    4,  40,   8,  30};
1252   Double_t minpho[9] = {  0.,  0., -0.1, 40,  0, -0.5,   0,-1.5,   0};
1253   Double_t maxpho[9] = {100., 30.,  0.5, 100,   1,  3.5,   2, 6.5,  30};
1254
1255   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);
1256   if(fqahist==1)fOutputList->Add(fInvmassLS);
1257
1258   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);
1259   if(fqahist==1)fOutputList->Add(fInvmassULS);
1260
1261   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);
1262   if(fqahist==1)fOutputList->Add(fInvmassLSmc);
1263
1264   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);
1265   if(fqahist==1)fOutputList->Add(fInvmassULSmc);
1266
1267   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 );
1268   fInvmassLSreco->Sumw2();
1269   fOutputList->Add(fInvmassLSreco);
1270
1271   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 );
1272   fInvmassULSreco->Sumw2();
1273   fOutputList->Add(fInvmassULSreco);
1274
1275   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 );
1276   fInvmassLSmc0->Sumw2();
1277   fOutputList->Add(fInvmassLSmc0);
1278
1279   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 );
1280   fInvmassLSmc1->Sumw2();
1281   fOutputList->Add(fInvmassLSmc1);
1282
1283   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 );
1284   fInvmassLSmc2->Sumw2();
1285   fOutputList->Add(fInvmassLSmc2);
1286
1287   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 );
1288   fInvmassLSmc3->Sumw2();
1289   fOutputList->Add(fInvmassLSmc3);
1290
1291   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 );
1292   fInvmassULSmc0->Sumw2();
1293   fOutputList->Add(fInvmassULSmc0);
1294
1295   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 );
1296   fInvmassULSmc1->Sumw2();
1297   fOutputList->Add(fInvmassULSmc1);
1298
1299   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 );
1300   fInvmassULSmc2->Sumw2();
1301   fOutputList->Add(fInvmassULSmc2);
1302
1303   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 );
1304   fInvmassULSmc3->Sumw2();
1305   fOutputList->Add(fInvmassULSmc3);
1306
1307   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1308   fOutputList->Add(fOpeningAngleLS);
1309
1310   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1311   fOutputList->Add(fOpeningAngleULS);
1312
1313   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
1314   fOutputList->Add(fPhotoElecPt);
1315
1316   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
1317   fOutputList->Add(fPhoElecPt);
1318
1319   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
1320   fOutputList->Add(fPhoElecPtM20);
1321
1322   fPhoElecPtM20Mass = new TH2F("fPhoElecPtM20Mass", "Pho-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1323   fPhoElecPtM20Mass->Sumw2();
1324   fOutputList->Add(fPhoElecPtM20Mass);
1325
1326   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
1327   fOutputList->Add(fSameElecPt);
1328
1329   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
1330   fOutputList->Add(fSameElecPtM20);
1331
1332   fSameElecPtM20Mass = new TH2F("fSameElecPtM20Mass", "Same-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1333   fSameElecPtM20Mass->Sumw2();
1334   fOutputList->Add(fSameElecPtM20Mass);
1335
1336   fSemiElecPtM20 = new TH2F("fSemiElecPtM20", "Semi-inclusive electron pt with M20",200,0,100,100,0,50);
1337   fOutputList->Add(fSemiElecPtM20);
1338
1339   fCent = new TH1F("fCent","Centrality",200,0,100) ;
1340   fOutputList->Add(fCent);
1341
1342   // Make common binning
1343   const Double_t kMinP = 0.;
1344   const Double_t kMaxP = 20.;
1345
1346   //+++ 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
1347   // 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)
1348   Int_t nBins[16] =  {  100,     7,  60,    20,    90,  100,   25,     60,    60,  100,    40,  10,  250,  100, 100,    8};
1349   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};
1350   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};
1351   //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);
1352   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);
1353   //if(fqahist==1)fOutputList->Add(fEleInfo);
1354
1355   // Make common binning
1356   Int_t nBinsEop[3] =  { 10, 50, 100};
1357   Double_t minEop[3] = {  0,  0,  -5};
1358   Double_t maxEop[3] = {100, 50,   5};
1359   fElenSigma= new THnSparseD("fElenSigma", "Electron nSigma; cent; pT [GeV/c]; nSigma", 3, nBinsEop, minEop, maxEop);
1360   fOutputList->Add(fElenSigma);
1361
1362
1363   //<---  Trigger info
1364   /*
1365   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
1366   fOutputList->Add(fClsEBftTrigCut);
1367
1368   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
1369   fOutputList->Add(fClsEAftTrigCut);
1370
1371   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
1372   fOutputList->Add(fClsEAftTrigCut1);
1373
1374   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
1375   fOutputList->Add(fClsEAftTrigCut2);
1376
1377   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
1378   fOutputList->Add(fClsEAftTrigCut3);
1379
1380   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
1381   fOutputList->Add(fClsEAftTrigCut4);
1382
1383   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1384   fOutputList->Add(fClsETime);
1385
1386   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1387   fOutputList->Add(fClsETime1);
1388
1389   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
1390   fOutputList->Add(fTrigTimes);
1391
1392   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
1393   fOutputList->Add(fCellCheck);
1394   */
1395   //<---------- MC
1396
1397   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1398   fOutputList->Add(fInputHFEMC);
1399
1400   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
1401   fOutputList->Add(fInputAlle);
1402
1403   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1404   fOutputList->Add(fIncpTMChfe);
1405
1406   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
1407   fOutputList->Add(fIncpTMChfeAll);
1408
1409   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1410   fOutputList->Add(fIncpTMCM20hfe);
1411
1412   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
1413   fOutputList->Add(fIncpTMCM20hfeAll);
1414
1415   fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1416   fOutputList->Add(fIncpTMCM20hfeCheck);
1417
1418   Int_t nBinspho2[5] =  { 200, 100,    7,    3, 700};
1419   Double_t minpho2[5] = {  0.,  0., -2.5, -0.5, 0.};
1420   Double_t maxpho2[5] = {100., 50.,  4.5,  2.5, 70.};
1421
1422   fInputHFEMC_weight = new THnSparseD("fInputHFEMC_weight", "MC HFE electron pt",5,nBinspho2,minpho2,maxpho2);
1423   fOutputList->Add(fInputHFEMC_weight);
1424
1425   fIncpTMCM20hfeCheck_weight = new THnSparseD("fIncpTMCM20hfeCheck_weight", "HFE electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1426   fOutputList->Add(fIncpTMCM20hfeCheck_weight);
1427
1428   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
1429   fOutputList->Add(fIncpTMCpho);
1430
1431   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1432   fOutputList->Add(fIncpTMCM20pho);
1433
1434   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1435   fOutputList->Add(fPhoElecPtMC);
1436
1437   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1438   fOutputList->Add(fPhoElecPtMCM20);
1439
1440   fPhoElecPtMCM20Mass = new TH2D("fPhoElecPtMCM20Mass", "MC Pho-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1441   fOutputList->Add(fPhoElecPtMCM20Mass);
1442
1443   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1444   fOutputList->Add(fSameElecPtMC);
1445
1446   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1447   fOutputList->Add(fSameElecPtMCM20);
1448
1449   fSameElecPtMCM20Mass = new TH2D("fSameElecPtMCM20Mass", "MC Same-inclusive electron pt with M20 Mass",200,0,100,100,0,50);
1450   fOutputList->Add(fSameElecPtMCM20Mass);
1451
1452   fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1453   fIncpTMCM20pho_pi0e->Sumw2();
1454   fOutputList->Add(fIncpTMCM20pho_pi0e);
1455
1456   fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1457   fPhoElecPtMCM20_pi0e->Sumw2();
1458   fOutputList->Add(fPhoElecPtMCM20_pi0e);
1459
1460   fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1461   fSameElecPtMCM20_pi0e->Sumw2();
1462   fOutputList->Add(fSameElecPtMCM20_pi0e);
1463  //
1464   fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1465   fIncpTMCM20pho_eta->Sumw2();
1466   fOutputList->Add(fIncpTMCM20pho_eta);
1467
1468   fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1469   fPhoElecPtMCM20_eta->Sumw2();
1470   fOutputList->Add(fPhoElecPtMCM20_eta);
1471
1472   fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1473   fSameElecPtMCM20_eta->Sumw2();
1474   fOutputList->Add(fSameElecPtMCM20_eta);
1475   // ------------
1476   fIncpTMCpho_pi0e_TPC = new THnSparseD("fIncpTMCpho_pi0e_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1477   fIncpTMCpho_pi0e_TPC->Sumw2();
1478   fOutputList->Add(fIncpTMCpho_pi0e_TPC);
1479
1480   fPhoElecPtMC_pi0e_TPC = new THnSparseD("fPhoElecPtMC_pi0e_TPC", "MC Pho-inclusive electron pt with  pi0->e",5,nBinspho2,minpho2,maxpho2);
1481   fPhoElecPtMC_pi0e_TPC->Sumw2();
1482   fOutputList->Add(fPhoElecPtMC_pi0e_TPC);
1483
1484   fSameElecPtMC_pi0e_TPC = new THnSparseD("fSameElecPtMC_pi0e_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1485   fSameElecPtMC_pi0e_TPC->Sumw2();
1486   fOutputList->Add(fSameElecPtMC_pi0e_TPC);
1487
1488   fIncpTMCpho_eta_TPC = new THnSparseD("fIncpTMCpho_eta_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1489   fIncpTMCpho_eta_TPC->Sumw2();
1490   fOutputList->Add(fIncpTMCpho_eta_TPC);
1491
1492   fPhoElecPtMC_eta_TPC = new THnSparseD("fPhoElecPtMC_eta_TPC", "MC Pho-inclusive electron pt with  pi0->e",5,nBinspho2,minpho2,maxpho2);
1493   fPhoElecPtMC_eta_TPC->Sumw2();
1494   fOutputList->Add(fPhoElecPtMC_eta_TPC);
1495
1496   fSameElecPtMC_eta_TPC = new THnSparseD("fSameElecPtMC_eta_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1497   fSameElecPtMC_eta_TPC->Sumw2();
1498   fOutputList->Add(fSameElecPtMC_eta_TPC);
1499
1500
1501   //-------------
1502
1503   CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1504   fOutputList->Add(CheckNclust);
1505
1506   CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1507   fOutputList->Add(CheckNits);
1508
1509   CheckDCA = new TH2D("CheckDCA","DCA check",200,-5,5,200,-5,5);
1510   fOutputList->Add(CheckDCA);
1511   /*
1512   Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1513   fOutputList->Add(Hpi0pTcheck);
1514
1515   HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1516   fOutputList->Add(HETApTcheck);
1517   */
1518
1519   Int_t nBinspho3[3] =  { 200, 100, 3};
1520   Double_t minpho3[3] = {  0.,  0., -0.5};
1521   Double_t maxpho3[3] = {100., 50., 2.5};
1522
1523   Hpi0pTcheck = new THnSparseD("Hpi0pTcheck","Pi0 pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1524   fOutputList->Add(Hpi0pTcheck);
1525
1526   HETApTcheck = new THnSparseD("HETApTcheck","Eta pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1527   fOutputList->Add(HETApTcheck);
1528   //--
1529   HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1530   fOutputList->Add(HphopTcheck);
1531   //
1532   HDpTcheck = new TH2D("HDpTcheck","D pT from Hijing",100,0,50,3,-0.5,2.5);
1533   fOutputList->Add(HDpTcheck);
1534
1535   HBpTcheck = new TH2D("HBpTcheck","B pT from Hijing",100,0,50,3,-0.5,2.5);
1536   fOutputList->Add(HBpTcheck);
1537   //
1538
1539   fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1540   fOutputList->Add(fpTCheck);
1541
1542   fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1543   fOutputList->Add(fMomDtoE);
1544
1545   fLabelCheck = new TH2D("fLabelCheck","MC label",50,0,50,5,-1.5,3.5);
1546   fOutputList->Add(fLabelCheck);
1547
1548   fgeoFake = new TH2D("fgeoFake","Label==0 eta and phi",628,0,6.28,200,-1,1);
1549   fOutputList->Add(fgeoFake);
1550
1551   fFakeTrk0 = new TH2D("fFakeTrk0","fake trakcs",10,0,100,20,0,20);
1552   fOutputList->Add(fFakeTrk0);
1553
1554   fFakeTrk1 = new TH2D("fFakeTrk1","true all e a.f. eID",10,0,100,20,0,20);
1555   fOutputList->Add(fFakeTrk1);
1556
1557   ftimingEle = new TH2D("ftimingEle","electron TOF",100,0,20,100,1e-7,1e-6);
1558   fOutputList->Add(ftimingEle);
1559
1560   // eta correction
1561   // note: parameters 01/31new.TPCnSigmaEtaDep
1562   // 70-90 delta_eta = 0.2
1563
1564   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};
1565   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)
1566   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)
1567   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)
1568   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)
1569   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)
1570   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)
1571   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)
1572
1573   fnSigEtaCorr[0] = new TGraphErrors(12,etaval,corr0); // 0-10
1574   fnSigEtaCorr[1] = new TGraphErrors(12,etaval,corr1); // 10-20
1575   fnSigEtaCorr[2] = new TGraphErrors(12,etaval,corr2); // 20-30
1576   fnSigEtaCorr[3] = new TGraphErrors(12,etaval,corr3); // 30-40
1577   fnSigEtaCorr[4] = new TGraphErrors(12,etaval,corr4); // 40-50
1578   fnSigEtaCorr[5] = new TGraphErrors(12,etaval,corr5); // 50-70
1579   fnSigEtaCorr[6] = new TGraphErrors(12,etaval,corr6); // 70-90
1580
1581   fIncMaxE = new TH2D("fIncMaxE","Inc",10,0,100,10,0,100);
1582   fOutputList->Add(fIncMaxE);
1583
1584   fIncReco = new TH2D("fIncReco","Inc",10,0,100,100,0,500);
1585   fOutputList->Add(fIncReco);
1586
1587   fPhoReco = new TH2D("fPhoReco","Pho",10,0,100,100,0,500);
1588   fOutputList->Add(fPhoReco);
1589
1590   fSamReco = new TH2D("fSamReco","Same",10,0,100,100,0,500);
1591   fOutputList->Add(fSamReco);
1592
1593   fIncRecoMaxE = new TH2D("fIncRecoMaxE","Inc",10,0,100,100,0,500);
1594   fOutputList->Add(fIncRecoMaxE);
1595
1596   fPhoRecoMaxE = new TH2D("fPhoRecoMaxE","Pho",10,0,100,100,0,500);
1597   fOutputList->Add(fPhoRecoMaxE);
1598
1599   fSamRecoMaxE = new TH2D("fSamRecoMaxE","Same",10,0,100,100,0,500);
1600   fOutputList->Add(fSamRecoMaxE);
1601
1602   fPhoVertexReco_TPC = new TH2D("fPhoVertexReco_TPC","photon production Vertex mass selection TPC",40,0,20,250,0,50);
1603   fPhoVertexReco_TPC->Sumw2();
1604   fOutputList->Add(fPhoVertexReco_TPC);
1605
1606   fPhoVertexReco_TPC_Invmass = new TH2D("fPhoVertexReco_TPC_Invmass","photon production Vertex mass selection TPC Invmass",40,0,20,250,0,50);
1607   fPhoVertexReco_TPC_Invmass->Sumw2();
1608   fOutputList->Add(fPhoVertexReco_TPC_Invmass);
1609
1610   fPhoVertexReco_EMCal = new TH2D("fPhoVertexReco_EMCal","photon production Vertex mass selection",40,0,20,250,0,50);
1611   fPhoVertexReco_EMCal->Sumw2();
1612   fOutputList->Add(fPhoVertexReco_EMCal);
1613
1614   fPhoVertexReco_Invmass = new TH2D("fPhoVertexReco_Invmass","photon production Vertex mass selection",40,0,20,250,0,50);
1615   fPhoVertexReco_Invmass->Sumw2();
1616   fOutputList->Add(fPhoVertexReco_Invmass);
1617
1618   fPhoVertexReco_step0= new TH2D("fPhoVertexReco_step0","photon production Vertex mass selection",40,0,20,250,0,50);
1619   fPhoVertexReco_step0->Sumw2();
1620   fOutputList->Add(fPhoVertexReco_step0);
1621
1622   fPhoVertexReco_step1= new TH2D("fPhoVertexReco_step1","photon production Vertex mass selection",40,0,20,250,0,50);
1623   fPhoVertexReco_step1->Sumw2();
1624   fOutputList->Add(fPhoVertexReco_step1);
1625
1626   fMatchV0_0 = new TH1D("fMatchV0_0","V0 match",100,0,20);
1627   fOutputList->Add(fMatchV0_0);
1628
1629   fMatchV0_1 = new TH1D("fMatchV0_1","V0 match",100,0,20);
1630   fOutputList->Add(fMatchV0_1);
1631
1632   fMatchMC_0 = new TH1D("fMatchMC_0","MC match",100,0,20);
1633   fOutputList->Add(fMatchMC_0);
1634
1635   fMatchMC_1 = new TH1D("fMatchMC_1","MC match",100,0,20);
1636   fOutputList->Add(fMatchMC_1);
1637
1638   fpair = new TH2D("fpair","pair of associate",100,0,20,21,-10.5,10.5);
1639   fOutputList->Add(fpair);
1640
1641   fFakeRejection0 = new TH2D("fFakeRejection0","TPC PID",2,-0.5,1.5,100,0,20);
1642   fFakeRejection0->Sumw2();
1643   fOutputList->Add(fFakeRejection0);
1644
1645   fFakeRejection1 = new TH2D("fFakeRejection1","TPC PID + Tr match",2,-0.5,1.5,100,0,20);
1646   fFakeRejection1->Sumw2();
1647   fOutputList->Add(fFakeRejection1);
1648
1649   fFakeRejection2 = new TH2D("fFakeRejection2","TPC PID + Tr match + E/p",2,-0.5,1.5,100,0,20);
1650   fFakeRejection2->Sumw2();
1651   fOutputList->Add(fFakeRejection2);
1652
1653   EopFake = new TH2D("EopFake","negative track Eop",20,0,20,200,0,2);
1654   EopFake->Sumw2();
1655   fOutputList->Add(EopFake);
1656
1657   EopTrue = new TH2D("EopTrue","true track Eop",20,0,20,200,0,2);
1658   EopTrue->Sumw2();
1659   fOutputList->Add(EopTrue);
1660
1661   MatchFake = new TH2D("MatchFake","negative track Match",20,0,20,100,0,0.05);
1662   MatchFake->Sumw2();
1663   fOutputList->Add(MatchFake);
1664
1665   MatchTrue = new TH2D("MatchTrue","true track Match",20,0,20,100,0,0.05);
1666   MatchTrue->Sumw2();
1667   fOutputList->Add(MatchTrue);
1668
1669   MatchTrCheck = new TH2D("MatchTrCheck","Check Tr Match",2010,-10,2000,2010,-10,2000);
1670   fOutputList->Add(MatchTrCheck);
1671
1672   MatchTrEop = new TH2D("MatchTrEop","Check Tr Match E/p",200,0,2,200,0,2);
1673   fOutputList->Add(MatchTrEop);
1674
1675   PostData(1,fOutputList);
1676 }
1677
1678 //________________________________________________________________________
1679 void AliAnalysisTaskHFECal::Terminate(Option_t *)
1680 {
1681   // Info("Terminate");
1682         AliAnalysisTaskSE::Terminate();
1683 }
1684
1685 //________________________________________________________________________
1686 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1687 {
1688   // Check single track cuts for a given cut step
1689   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1690   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1691   return kTRUE;
1692 }
1693 //_________________________________________
1694 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)
1695 {
1696   //Identify non-heavy flavour electrons using Invariant mass method
1697
1698   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1699   fTrackCuts->SetRequireTPCRefit(kTRUE);
1700   fTrackCuts->SetRequireITSRefit(kTRUE);
1701   fTrackCuts->SetEtaRange(-0.9,0.9);
1702   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1703   fTrackCuts->SetMinNClustersTPC(90);
1704
1705   //const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1706   Double_t bfield = fESD->GetMagneticField();
1707   Double_t emass = 0.51*0.001; // (0.51 MeV)
1708
1709   double ptEle = track->Pt();  //add
1710   if(ibgevent==0 && w > 0.0)
1711      {
1712       fpTCheck->Fill(ptEle,w);
1713      }
1714
1715   Bool_t flagPhotonicElec = kFALSE;
1716   Bool_t flagConvinatElec = kFALSE;
1717
1718   int p1 = 0;
1719   if(mce==3)
1720      {
1721        Int_t label = TMath::Abs(track->GetLabel());
1722        TParticle* particle = stack->Particle(label);
1723        p1 = particle->GetFirstMother();
1724      }
1725
1726   int numULS = 0;
1727   int numLS = 0;
1728
1729   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1730     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1731     if (!trackAsso) {
1732       printf("ERROR: Could not receive track %d\n", jTracks);
1733       continue;
1734     }
1735     if(itrack==jTracks)continue;
1736     int jbgevent = 0;
1737
1738     int p2 = 0;
1739     if(mce==3)
1740     {
1741       Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1742       TParticle* particle2 = stack->Particle(label2);
1743       Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1744       if(MChijing_ass)jbgevent =1;
1745       if(particle2->GetFirstMother()>-1)
1746          p2 = particle2->GetFirstMother();
1747     }
1748
1749     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1750     Double_t mass=999., width = -999;
1751     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1752
1753     //ptPrim = track->Pt();
1754     ptPrim = ptEle;
1755
1756     dEdxAsso = trackAsso->GetTPCsignal();
1757     ptAsso = trackAsso->Pt();
1758     Int_t chargeAsso = trackAsso->Charge();
1759     Int_t charge = track->Charge();
1760
1761
1762     if(ptAsso <fMimpTassCut) continue;
1763     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1764     double fTPCnSigmaAss = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
1765     //cout << "fTPCnSigmaAss = " << fTPCnSigmaAss << endl;
1766     //cout << "fTPCnSigmaAss Cut = " << fMimNsigassCut << endl;
1767     if(fTPCnSigmaAss <fMimNsigassCut || fTPCnSigmaAss>5) continue;
1768     //cout << "fTPCnSigmaAss a.f. cut = " << fTPCnSigmaAss << endl;
1769
1770     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1771     if(charge>0) fPDGe1 = -11;
1772     if(chargeAsso>0) fPDGe2 = -11;
1773
1774     //printf("chargeAsso = %d\n",chargeAsso);
1775     //printf("charge = %d\n",charge);
1776     if(charge == chargeAsso) fFlagLS = kTRUE;
1777     if(charge != chargeAsso) fFlagULS = kTRUE;
1778
1779     //printf("fFlagLS = %d\n",fFlagLS);
1780     //printf("fFlagULS = %d\n",fFlagULS);
1781     //printf("\n");
1782
1783     AliKFParticle::SetField(bfield);
1784     AliKFParticle ge1(*track, fPDGe1);
1785     AliKFParticle ge2(*trackAsso, fPDGe2);
1786
1787     if(fSetMassNonlinear)
1788       {
1789        ge1.SetNonlinearMassConstraint(emass);
1790        ge2.SetNonlinearMassConstraint(emass);
1791       }
1792
1793     AliKFParticle recg(ge1, ge2);
1794
1795     /*
1796     // vertex
1797     AliKFVertex primV(*pVtx);
1798     primV += recg;
1799     primV -= ge1;
1800     primV -= ge2;
1801     recg.SetProductionVertex(primV);
1802     */
1803
1804     // check chi2
1805     if(recg.GetNDF()<1) continue;
1806     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1807     Double_t chi2cut = 3.0;
1808
1809     // mass.
1810     if(fSetMassConstraint)
1811       {
1812        recg.SetMassConstraint(0,0.0001);
1813        chi2cut = 30.0;
1814       }
1815     recg.GetMass(mass,width);
1816
1817     if(fSetMassWidthCut && width>1e10)continue;
1818
1819         // angle
1820     openingAngle = ge1.GetAngle(ge2);
1821     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1822     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1823
1824     double ishower = 0;
1825     if(shower>0.0 && shower<0.3)ishower = 1;
1826
1827     if(iCal==1)
1828       {
1829         double phoinfo[9];
1830         phoinfo[0] = cent;
1831         phoinfo[1] = ptPrim;
1832         phoinfo[2] = mass;
1833         phoinfo[3] = nSig;
1834         phoinfo[4] = openingAngle;
1835         phoinfo[5] = ishower;
1836         phoinfo[6] = ep;
1837         phoinfo[7] = mce;
1838         phoinfo[8] = ptAsso;
1839
1840         if(fFlagLS) fInvmassLS->Fill(phoinfo);
1841         if(fFlagULS) fInvmassULS->Fill(phoinfo);
1842         if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1843         if(fFlagULS && ibgevent==0 && jbgevent==0)
1844         {
1845               fInvmassULSmc->Fill(phoinfo,w);
1846         }
1847        //printf("fInvmassCut %f\n",fInvmassCut);
1848        //printf("openingAngle %f\n",fOpeningAngleCut);
1849       }
1850
1851     // angle cut
1852     if(openingAngle > fOpeningAngleCut) continue;
1853     // chi2 cut
1854     //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue;
1855     if(chi2recg>chi2cut) continue;
1856
1857     if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass);
1858     if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass);
1859
1860     // check double count
1861     if(mass<fInvmassCut && fFlagULS)numULS++;
1862     if(mass<fInvmassCut && fFlagLS)numLS++;
1863
1864     // for real data
1865     //printf("mce =%f\n",mce);
1866
1867     if(mce<-0.5) // mce==-1. is real
1868     {
1869             //printf("Real data\n");
1870
1871             // count
1872             if(mass<fInvmassCut && fFlagULS && shower>0.0 && shower<0.3 && iCal==1)fPhoElecPtM20Mass->Fill(cent,ptPrim);
1873             if(mass<fInvmassCut && fFlagLS && shower>0.0 && shower<0.3 && iCal==1)fSameElecPtM20Mass->Fill(cent,ptPrim);
1874
1875             // flag
1876             if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1877                     flagPhotonicElec = kTRUE;
1878             }
1879             if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1880                     flagConvinatElec = kTRUE;
1881             }
1882     }
1883     // for MC data
1884     else
1885     {
1886             // count
1887             if(mass<fInvmassCut && fFlagULS && shower>0.0 && shower<0.3 && mce>2.1 && iCal==1)fPhoElecPtMCM20Mass->Fill(cent,ptPrim,w);
1888             if(mass<fInvmassCut && fFlagLS && shower>0.0 && shower<0.3 && mce>2.1 && iCal==1)fSameElecPtMCM20Mass->Fill(cent,ptPrim,w);
1889
1890             //printf("MC data\n");
1891
1892             if(w>0.0)
1893             {
1894                     //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1895                     if(iCal==1)
1896                       {
1897                        if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass);
1898                        if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass);
1899                        if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass);
1900                        if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass);
1901                        if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass);
1902                        if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass);
1903                        if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass);
1904                        if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass);
1905                       }
1906             }
1907
1908             if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1909                     //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ //<--- only MC train (55,56) v5-03-68-AN , 69 & v5-05-70-AN (till 74)
1910                     flagPhotonicElec = kTRUE;
1911             }
1912             if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1913                     flagConvinatElec = kTRUE;
1914             }
1915     }
1916
1917   } // end of associate loop
1918
1919   if(numULS>0 || numLS>0)
1920     {
1921      int numPair = numULS-numLS;
1922      if(iCal==1)fpair->Fill(ptEle,numPair);
1923     }
1924
1925   fFlagPhotonic = flagPhotonicElec;
1926   fFlagConvinat = flagConvinatElec;
1927
1928 }
1929
1930 //-------------------------------------------
1931
1932 void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1933 {
1934  //int label = 99999;
1935  //int pid = 99999;
1936
1937  if(part->GetFirstMother()>-1)
1938    {
1939     label = part->GetFirstMother();
1940     pid = stack->Particle(label)->GetPdgCode();
1941    }
1942    //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
1943 }
1944
1945 double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1946 {
1947         double weight = 1.0;
1948
1949         if(mcPi0pT>0.0 && mcPi0pT<5.0)
1950         {
1951                 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1952         }
1953         else
1954         {
1955                 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
1956         }
1957   return weight;
1958 }
1959
1960 double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
1961 {
1962   double weight = 1.0;
1963
1964   weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
1965   return weight;
1966 }
1967
1968
1969 //_________________________________________
1970 void AliAnalysisTaskHFECal::FindTriggerClusters()
1971 {
1972   //cout << "finding trigger patch" << endl;
1973   // constants
1974   const int nModuleCols = 2;
1975   const int nModuleRows = 5;
1976   const int nColsFeeModule = 48;
1977   const int nRowsFeeModule = 24;
1978   const int nColsFaltroModule = 24;
1979   const int nRowsFaltroModule = 12;
1980   //const int faltroWidthMax = 20;
1981
1982   // part 1, trigger extraction -------------------------------------
1983   Int_t globCol, globRow;
1984   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
1985   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
1986
1987   //Int_t trigtimes[faltroWidthMax];
1988   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1989   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
1990   //Double_t fTrigCutLow = 6;
1991   //Double_t fTrigCutHigh = 10;
1992   Double_t fTimeCutLow =  469e-09;
1993   Double_t fTimeCutHigh = 715e-09;
1994
1995   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
1996
1997   // erase trigger maps
1998   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
1999   {
2000     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
2001     {
2002       ftriggersCut[i][j] = 0;
2003       ftriggers[i][j] = 0;
2004       ftriggersTime[i][j] = 0;
2005     }
2006   }
2007
2008   Int_t iglobCol=0, iglobRow=0;
2009   // go through triggers
2010   if( fCaloTrigger->GetEntries() > 0 )
2011   {
2012     // needs reset
2013     fCaloTrigger->Reset();
2014     while( fCaloTrigger->Next() )
2015     {
2016       fCaloTrigger->GetPosition( globCol, globRow );
2017       fCaloTrigger->GetNL0Times( ntimes );
2018       /*
2019       // no L0s
2020       if( ntimes < 1 )   continue;
2021       // get precise timings
2022       fCaloTrigger->GetL0Times( trigtimes );
2023       trigInCut = 0;
2024       for(Int_t i = 0; i < ntimes; i++ )
2025       {
2026         // save the first trigger time in channel
2027         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
2028           triggersTime[globCol][globRow] = trigtimes[i];
2029         //printf("trigger times: %d\n",trigtimes[i]);
2030         // check if in cut
2031         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
2032           trigInCut = 1;
2033
2034         fTrigTimes->Fill(trigtimes[i]);
2035       }
2036       */
2037
2038       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
2039       Int_t bit = 0;
2040       fCaloTrigger->GetTriggerBits(bit);
2041       //cout << "bit = " << bit << endl;
2042
2043       Int_t ts = 0;
2044       fCaloTrigger->GetL1TimeSum(ts);
2045       //cout << "ts = " << ts << endl;
2046       if (ts > 0)ftriggers[globCol][globRow] = 1;
2047       // number of triggered channels in event
2048       nTrigChannel++;
2049       // ... inside cut
2050       if(ts>0 && (bit >> 6 & 0x1))
2051       {
2052         iglobCol = globCol;
2053         iglobRow = globRow;
2054         nTrigChannelCut++;
2055         //cout << "ts cut = " << ts << endl;
2056         //cout << "globCol = " << globCol << endl;
2057         //cout << "globRow = " << globRow << endl;
2058         ftriggersCut[globCol][globRow] = 1;
2059       }
2060
2061     } // calo trigger entries
2062   } // has calo trigger entries
2063
2064   // part 2 go through the clusters here -----------------------------------
2065   //cout << " part 2 go through the clusters here ----------------------------------- " << endl;
2066   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
2067   Short_t cellAddr, nSACell;
2068   Int_t mclabel;
2069   //Int_t nSACell, iSACell, mclabel;
2070   Int_t iSACell;
2071   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
2072   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
2073
2074   TRefArray *fCaloClusters = new TRefArray();
2075   fESD->GetEMCALClusters( fCaloClusters );
2076   nCluster = fCaloClusters->GetEntries();
2077
2078
2079   // save all cells times if there are clusters
2080   if( nCluster > 0 ){
2081     // erase time map
2082     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
2083       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
2084         cellTime[i][j] = 0.;
2085         cellEnergy[i][j] = 0.;
2086       }
2087     }
2088
2089     // get cells
2090     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
2091     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
2092     nSACell = fCaloCells->GetNumberOfCells();
2093     for(iSACell = 0; iSACell < nSACell; iSACell++ ){
2094       // get the cell info *fCal
2095       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
2096
2097       // get cell position
2098       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
2099       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2100
2101       // convert co global phi eta
2102       gphi = iphi + nRowsFeeModule*(nSupMod/2);
2103       geta = ieta + nColsFeeModule*(nSupMod%2);
2104
2105       // save cell time and energy
2106       cellTime[geta][gphi] = cellTimeT;
2107       cellEnergy[geta][gphi] = cellAmp;
2108
2109     }
2110   }
2111
2112   Int_t nClusterTrig, nClusterTrigCut;
2113   UShort_t *cellAddrs;
2114   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
2115   Float_t clsPos[3] = {0.,0.,0.};
2116
2117   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
2118   {
2119     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2120     if(!cluster || !cluster->IsEMCAL()) continue;
2121
2122     // get cluster cells
2123     nCell = cluster->GetNCells();
2124
2125     // get cluster energy
2126     clsE = cluster->E();
2127
2128     // get cluster position
2129     cluster->GetPosition(clsPos);
2130     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2131     clsEta = clsPosVec.Eta();
2132     clsPhi = clsPosVec.Phi();
2133
2134     // get the cell addresses
2135     cellAddrs = cluster->GetCellsAbsId();
2136
2137     // check if the cluster contains cell, that was marked as triggered
2138     nClusterTrig = 0;
2139     nClusterTrigCut = 0;
2140
2141     // loop the cells to check, if cluser in acceptance
2142     // any cluster with a cell outside acceptance is not considered
2143     for( iCell = 0; iCell < nCell; iCell++ )
2144     {
2145      // check hot cell
2146      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
2147
2148       // get cell position
2149       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
2150       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2151
2152       // convert co global phi eta
2153       gphi = iphi + nRowsFeeModule*(nSupMod/2);
2154       geta = ieta + nColsFeeModule*(nSupMod%2);
2155
2156       if( cellTime[geta][gphi] > 0. ){
2157         clusterTime += cellTime[geta][gphi];
2158         gCell++;
2159       }
2160
2161       // get corresponding FALTRO
2162       fphi = gphi / 2;
2163       feta = geta / 2;
2164
2165         //cout << "fphi = " << fphi << endl;
2166         //cout << "feta = " << feta << endl;
2167
2168       // try to match with a triggered
2169       if( ftriggers[feta][fphi]==1)
2170       {  nClusterTrig++;
2171       }
2172       if( ftriggersCut[feta][fphi]==1)
2173       { nClusterTrigCut++;
2174       }
2175
2176       //cout << "nClusterTrigCut : " << nClusterTrigCut << endl;
2177
2178     } // cells
2179
2180
2181     if( gCell > 0 )
2182       clusterTime = clusterTime / (Double_t)gCell;
2183     // fix the reconstriction code time 100ns jumps
2184     if( fESD->GetBunchCrossNumber() % 4 < 2 )
2185       clusterTime -= 0.0000001;
2186
2187     //fClsETime->Fill(clsE,clusterTime);
2188     //fClsEBftTrigCut->Fill(clsE);
2189
2190     if(nClusterTrig>0){
2191       //fClsETime1->Fill(clsE,clusterTime);
2192     }
2193
2194     if(nClusterTrig>0){
2195       cluster->SetChi2(1);
2196       //fClsEAftTrigCut1->Fill(clsE);
2197     }
2198
2199     if(nClusterTrigCut>0){
2200       cluster->SetChi2(2);
2201       //fClsEAftTrigCut2->Fill(clsE);
2202     }
2203
2204     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
2205     {
2206       cluster->SetChi2(3);
2207       //fClsEAftTrigCut3->Fill(clsE);
2208     }
2209
2210     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
2211     {
2212       // cluster->SetChi2(4);
2213       //fClsEAftTrigCut4->Fill(clsE);
2214     }
2215     if(nClusterTrigCut<1)
2216     {
2217       cluster->SetChi2(0);
2218
2219       //fClsEAftTrigCut->Fill(clsE);
2220     }
2221
2222   } // clusters
2223 }
2224
2225 // <-------- only MC correction
2226 double AliAnalysisTaskHFECal::MCEopMeanCorrection(double pTmc, float central)
2227 {
2228   TF1 *fcorr0 = new TF1("fcorr0","[0]*tanh([1]+[2]*x)");
2229   TF1 *fcorr1 = new TF1("fcorr1","[0]*tanh([1]+[2]*x)");
2230
2231  double shift = 0.0;
2232
2233  if(central>0 && central<=10)
2234    {
2235     fcorr0->SetParameters(1.045,1.288,3.18e-01); //
2236     fcorr1->SetParameters(9.91e-01,3.466,2.344);
2237    }
2238  else if(central>10 && central<=20)
2239    {
2240     fcorr0->SetParameters(1.029,8.254e-01,4.07e-01);
2241     fcorr1->SetParameters(0.975,2.276,1.501e-01);
2242    }
2243  else if(central>20 && central<=30)
2244    {
2245     fcorr0->SetParameters(1.01,8.795e-01,3.904e-01);
2246     fcorr1->SetParameters(9.675e-01,1.654,2.583e-01);
2247    }
2248  else if(central>30 && central<=40)
2249    {
2250     fcorr0->SetParameters(1.00,1.466,2.305e-1);
2251     fcorr1->SetParameters(9.637e-01,1.473,2.754e-01);
2252    }
2253  else if(central>40 && central<=50)
2254    {
2255     fcorr0->SetParameters(1.00,1.422,1.518e-01);
2256     fcorr1->SetParameters(9.59e-01,1.421,2.931e-01);
2257    }
2258
2259  else if(central>50 && central<=70)
2260    {
2261     fcorr0->SetParameters(0.989,2.495,2.167);
2262     fcorr1->SetParameters(0.961,1.734,1.438e-01);
2263    }
2264  else if(central>70 && central<=100)
2265    {
2266     fcorr0->SetParameters(0.981,-3.373,3.93327);
2267     fcorr1->SetParameters(9.574e-01,1.698,1.58e-01);
2268    }
2269
2270
2271  shift = fcorr0->Eval(pTmc)-fcorr1->Eval(pTmc);
2272
2273  fcorr0->Delete();
2274  fcorr1->Delete();
2275
2276  return shift;
2277 }
2278
2279 // <-------- only Data correction
2280 double AliAnalysisTaskHFECal::NsigmaCorrection(double tmpeta, float central)
2281 {
2282  int icent = 0;
2283
2284  if(central>=0 && central<10)
2285    {
2286     icent = 0;
2287    }
2288  else if(central>=10 && central<20)
2289   {
2290    icent = 1;
2291   }
2292  else if(central>=20 && central<30)
2293   {
2294    icent = 2;
2295   }
2296  else if(central>=30 && central<40)
2297   {
2298    icent = 3;
2299   }
2300  else if(central>=40 && central<50)
2301   {
2302    icent = 4;
2303   }
2304  else if(central>=50 && central<70)
2305   {
2306    icent = 5;
2307   }
2308  else
2309   {
2310    icent = 6;
2311   }
2312
2313  double shift = fnSigEtaCorr[icent]->Eval(tmpeta);
2314
2315  //cout << "eta correction"<< endl;
2316  //cout << "cent = "<< central<< endl;
2317  //cout << "icent = "<< icent << endl;
2318  //cout << "shift = "<< shift << endl;
2319
2320  return shift;
2321
2322 }
2323