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