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