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