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