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