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
 714     int nITS = track->GetNcls(0);
 715
 716     fTrackPtBefTrkCuts->Fill(track->Pt());
 717
 718     UChar_t itsPixel = track->GetITSClusterMap();
 719     /*
 720     cout << "nITS = " << nITS << endl;
 721     if(itsPixel & BIT(0))cout << "1st layer hit" << endl;
 722     if(itsPixel & BIT(1))cout << "2nd layer hit" << endl;
 723     if(itsPixel & BIT(2))cout << "3rd layer hit" << endl;
 724     if(itsPixel & BIT(3))cout << "4th layer hit" << endl;
 725     if(itsPixel & BIT(4))cout << "5th layer hit" << endl;
 726     */
 727
 728     Bool_t kAnyITS = kFALSE;
 729     if((itsPixel & BIT(0)) || (itsPixel & BIT(1)))kAnyITS = kTRUE;
 730
 731
 732     if(mcPho)fPhoVertexReco_step0->Fill(track->Pt(),conv_proR); // check MC vertex
 733
 734     // RecKine: ITSTPC cuts
 735     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 736     if(mcPho)fPhoVertexReco_step1->Fill(track->Pt(),conv_proR); // check MC vertex
 737
 738     //RecKink
 739     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 740       if(track->GetKinkIndex(0) != 0) continue;
 741     }
 742     if(mcPho)fPhoVertexReco_step2->Fill(track->Pt(),conv_proR); // check MC vertex
 743
 744     // RecPrim
 745     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 746     if(mcPho)fPhoVertexReco_step3->Fill(track->Pt(),conv_proR); // check MC vertex
 747
 748     // HFEcuts: ITS layers cuts
 749     //if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 750     if(!kAnyITS) continue; // select ITS condition without HFE
 751     if(mcPho)fPhoVertexReco_step4->Fill(track->Pt(),conv_proR); // check MC vertex
 752
 753     // HFE cuts: TPC PID cleanup
 754     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 755     if(mcPho)fPhoVertexReco_step5->Fill(track->Pt(),conv_proR); // check MC vertex
 756
 757     int nTPCcl = track->GetTPCNcls();
 758     //int nTPCclF = track->GetTPCNclsF(); // warnings
 759     //int nITS = track->GetNcls(0);
 760
 761     if(mcPho)fPhoVertexReco_HFE->Fill(track->Pt(),conv_proR); // check MC vertex
 762
 763     fTrackPtAftTrkCuts->Fill(track->Pt());
 764
 765     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
 766     pt = track->Pt();
 767     if(pt<2.5)continue;
 768
 769     //Int_t charge = track->Charge();
 770     fTrkpt->Fill(pt);
 771     mom = track->P();
 772     phi = track->Phi();
 773     eta = track->Eta();
 774     float dca_xy;
 775     float dca_z;
 776     track->GetImpactParameters(dca_xy,dca_z);
 777
 778     dEdx = track->GetTPCsignal();
 779     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 780
 781     //cout << "nSigma correctoon-----" << endl;
 782     //cout << "org = " << fTPCnSigma << endl;
 783     /*
 784     if(!fmcData) // nsigma eta correction
 785        {
 786         double nSigexpCorr = NsigmaCorrection(eta,cent);
 787         fTPCnSigma -= nSigexpCorr;
 788        }
 789     */
 790     //cout << "correction = " << fTPCnSigma << endl;
 791
 792     //--- track cluster match
 793
 794         double ncells = -1.0;
 795         double m20 = -1.0;
 796         double m02 = -1.0;
 797         double disp = -1.0;
 798         double rmatch = -1.0;
 799         double nmatch = -1.0;
 800         //double oppstatus = 0.0;
 801         double emctof = 0.0;
 802         Bool_t MaxEmatch = kFALSE;
 803
 804     Int_t clsId = track->GetEMCALcluster();
 805     if (clsId>0){
 806       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
 807       if(clust && clust->IsEMCAL()){
 808
 809                 double clustE = clust->E();
 810                 if(clustE==maxE)MaxEmatch = kTRUE;
 811                 eop = clustE/fabs(mom);
 812
 813                 //double clustT = clust->GetTOF();
 814                 ncells = clust->GetNCells();
 815                 m02 = clust->GetM02();
 816                 m20 = clust->GetM20();
 817                 disp = clust->GetDispersion();
 818                 double delphi = clust->GetTrackDx();
 819                 double deleta = clust->GetTrackDz();
 820                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
 821                 nmatch = clust->GetNTracksMatched();
 822                 emctof = clust->GetTOF();
 823                 //cout << "emctof = " << emctof << endl;
 824                 cout << "eop org = "<< eop << endl;
 825                 double eoporg =  eop;
 826                 if(fmcData)
 827                   {
 828                    double mceopcorr = MCEopMeanCorrection(pt,cent);
 829                    eop += mceopcorr;
 830                   }
 831                 cout << "eop corr = " << eop << endl;
 832
 833                   double valdedx[16];
 834                   valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma;
 835                   //valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells,  valdedx[8] = nmatch; valdedx[9] = m20; valdedx[10] = mcpT;
 836                   valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = dca_xy,  valdedx[8] = dca_z; valdedx[9] = m20; valdedx[10] = mcpT;
 837                   valdedx[11] = cent; valdedx[12] = dEdx; valdedx[13] = eoporg; valdedx[14] = nTPCcl;
 838                   valdedx[15] = mcele;
 839                   fEleInfo->Fill(valdedx);
 840
 841       }
 842     }
 843
 844     //Get Cal info PID response
 845     /*
 846     double eop2;
 847     double ss[4];
 848     Double_t nSigmaEop = fPID->GetPIDResponse()->NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss);
 849     if(fTPCnSigma>-1.5 && fTPCnSigma<3.0 && nITS>2.5 && nTPCcl>100)
 850       {
 851        double valEop[3];
 852        valEop[0] = cent;
 853        valEop[1] = pt;
 854        valEop[2] = nSigmaEop;
 855        fElenSigma->Fill(valEop);
 856       }
 857     */
 858    // --- tarck cut & e ID
 859
 860     if(nITS<2.5)continue;
 861     if(nTPCcl<100)continue;
 862
 863     CheckNclust->Fill(nTPCcl);
 864     CheckNits->Fill(nITS);
 865     CheckDCA->Fill(dca_xy,dca_z);
 866     // check production vertex of photons
 867
 868
 869     fdEdxBef->Fill(mom,fTPCnSigma);
 870     fTPCnsigma->Fill(mom,fTPCnSigma);
 871     if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
 872
 873     //--- track accepted by HFE
 874
 875     Int_t pidpassed = 1;
 876     AliHFEpidObject hfetrack;
 877     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 878     hfetrack.SetRecTrack(track);
 879     double binct = 10.5;
 880     if((0.0< cent) && (cent<5.0)) binct = 0.5;
 881     if((5.0< cent) && (cent<10.0)) binct = 1.5;
 882     if((10.0< cent) && (cent<20.0)) binct = 2.5;
 883     if((20.0< cent) && (cent<30.0)) binct = 3.5;
 884     if((30.0< cent) && (cent<40.0)) binct = 4.5;
 885     if((40.0< cent) && (cent<50.0)) binct = 5.5;
 886     if((50.0< cent) && (cent<60.0)) binct = 6.5;
 887     if((60.0< cent) && (cent<70.0)) binct = 7.5;
 888     if((70.0< cent) && (cent<80.0)) binct = 8.5;
 889     if((80.0< cent) && (cent<90.0)) binct = 9.5;
 890     if((90.0< cent) && (cent<100.0)) binct = 10.5;
 891
 892     hfetrack.SetCentrality((int)binct); //added
 893     hfetrack.SetPbPb();
 894     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 895
 896     if(pidpassed==0) continue;
 897
 898     //--- photonic ID
 899
 900     Bool_t fFlagPhotonicElec = kFALSE;
 901     Bool_t fFlagConvinatElec = kFALSE;
 902
 903     if(fSetKFpart)
 904       {
 905        SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta);
 906       }
 907     else
 908      {
 909       SelectPhotonicElectron2(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta);
 910       }
 911
 912     //--- check reco eff. with TPC
 913     double phoval[5];
 914     phoval[0] = cent;
 915     phoval[1] = pt;
 916     phoval[2] = fTPCnSigma;
 917     phoval[3] = iHijing;
 918     phoval[4] = mcMompT;
 919
 920     //if((fTPCnSigma >= -5.0 && fTPCnSigma <= 5) && mcele>-1 && mcPho && mcOrgPi0)
 921     if((fTPCnSigma >= -5.0 && fTPCnSigma <= 5) && (mcOrgPi0 || mcOrgEta))
 922       {
 923         if(iHijing==1)mcWeight = 1.0;
 924         if(mcOrgPi0)
 925           {
 926            fIncpTMCpho_pi0e_TPC->Fill(phoval,mcWeight);
 927            if(fFlagPhotonicElec) fPhoElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
 928            if(fFlagConvinatElec) fSameElecPtMC_pi0e_TPC->Fill(phoval,mcWeight);
 929           }
 930         if(mcOrgEta)
 931           {
 932            fIncpTMCpho_eta_TPC->Fill(phoval,mcWeight);
 933            if(fFlagPhotonicElec) fPhoElecPtMC_eta_TPC->Fill(phoval,mcWeight);
 934            if(fFlagConvinatElec) fSameElecPtMC_eta_TPC->Fill(phoval,mcWeight);
 935           }
 936       }
 937
 938     //+++++++  E/p cut ++++++++++++++++
 939
 940     if(eop<0.9 || eop>1.3)continue;
 941
 942     fTrkEovPAft->Fill(pt,eop);
 943     fdEdxAft->Fill(mom,fTPCnSigma);
 944
 945     // Fill real data
 946     fIncpT->Fill(cent,pt);
 947     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
 948     if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt);
 949
 950     if(m20>0.0 && m20<0.3)
 951       {
 952        fIncpTM20->Fill(cent,pt);
 953        ftimingEle->Fill(pt,emctof);
 954        if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt);
 955        if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt);
 956        if(!fFlagPhotonicElec) fSemiElecPtM20->Fill(cent,pt);
 957      }
 958
 959
 960     //--------
 961     /*
 962     double recopT =  SumpT(iTracks,track);
 963
 964     if(m20>0.0 && m20<0.3)
 965       {
 966        if(MaxEmatch)fIncMaxE->Fill(cent,pt);
 967        if(pt>5.0)
 968          {
 969           fIncReco->Fill(cent,recopT);
 970           if(fFlagPhotonicElec) fPhoReco->Fill(cent,recopT);
 971           if(fFlagConvinatElec) fSamReco->Fill(cent,recopT);
 972           if(MaxEmatch)
 973             {
 974              fIncRecoMaxE->Fill(cent,recopT);
 975              if(fFlagPhotonicElec) fPhoRecoMaxE->Fill(cent,recopT);
 976              if(fFlagConvinatElec) fSamRecoMaxE->Fill(cent,recopT);
 977             }
 978          }
 979      }
 980     */
 981
 982     // MC
 983     // check label for electron candidiates
 984
 985     int idlabel = 1;
 986     if(mcLabel==0)idlabel = 0;
 987     fLabelCheck->Fill(pt,idlabel);
 988     if(mcLabel==0)fgeoFake->Fill(phi,eta);
 989
 990     if(mcLabel<0 && m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)
 991        {
 992         fFakeTrk0->Fill(cent,pt);
 993        }
 994
 995     if(mcele>-1) // select MC electrons
 996       {
 997
 998           fIncpTMChfeAll->Fill(cent,pt);
 999           if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt);
1000           if(m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)fFakeTrk1->Fill(cent,pt);
1001
1002        if(mcBtoE || mcDtoE) // select B->e & D->e
1003          {
1004           fIncpTMChfe->Fill(cent,pt);
1005           if(m20>0.0 && m20<0.3)
1006             {
1007                 //cout << "MC label = " << mcLabel << endl;
1008                 fIncpTMCM20hfe->Fill(cent,pt);
1009                 fIncpTMCM20hfeCheck->Fill(cent,mcpT);
1010                 fIncpTMCM20hfeCheck_weight->Fill(phoval);
1011             }
1012          }
1013
1014        if(mcPho) // select photonic electrons
1015         {
1016
1017          fIncpTMCpho->Fill(phoval);
1018          if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval);
1019          if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval);
1020
1021          if(m20>0.0 && m20<0.3)
1022            {
1023             fIncpTMCM20pho->Fill(phoval);
1024             if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval);
1025             if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval);
1026             // pi0->g->e
1027             if(mcWeight>-1)
1028               {
1029                if(iHijing==1)mcWeight = 1.0;
1030                if(mcOrgPi0)
1031                  {
1032                   fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight);
1033                   if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight);
1034                   if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight);
1035
1036                   // check production vertex
1037                   fPhoVertexReco_EMCal->Fill(track->Pt(),conv_proR);
1038                   if(fFlagPhotonicElec)fPhoVertexReco_Invmass->Fill(track->Pt(),conv_proR);
1039
1040                   }
1041                // --- eta
1042                if(mcOrgEta)
1043                  {
1044                   fIncpTMCM20pho_eta->Fill(phoval,mcWeight);
1045                   if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight);
1046                   if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight);
1047                  }
1048                 // check production vertex
1049               }
1050            }
1051         }
1052       }
1053  }
1054  PostData(1, fOutputList);
1055 }
1056 //_________________________________________
1057 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
1058 {
1059   //--- Check MC
1060
1061   //Bool_t mcData = kFALSE;
1062   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
1063     {
1064      fmcData = kTRUE;
1065      printf("+++++ MC Data available");
1066     }
1067   if(fmcData)
1068     {
1069      printf("++++++++= MC analysis \n");
1070     }
1071   else
1072    {
1073      printf("++++++++ real data analysis \n");
1074    }
1075
1076    printf("+++++++ QA hist %d",fqahist);
1077
1078   //---- Geometry
1079   fGeom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
1080
1081   //--------Initialize PID
1082   //fPID->SetHasMCData(kFALSE);
1083   fPID->SetHasMCData(fmcData);
1084   if(!fPID->GetNumberOfPIDdetectors())
1085     {
1086       fPID->AddDetector("TPC", 0);
1087       //fPID->AddDetector("EMCAL", 1); <--- apply PID selection in task
1088     }
1089
1090   Double_t params[4];
1091   const char *cutmodel;
1092   cutmodel = "pol0";
1093   params[0] = -1.0; //sigma min
1094   double maxnSig = 3.0;
1095   if(fmcData)
1096     {
1097      params[0] = -5.0; //sigma min
1098      maxnSig = 5.0;
1099     }
1100
1101   for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig);
1102
1103   fPID->SortDetectors();
1104   fPIDqa = new AliHFEpidQAmanager();
1105   fPIDqa->Initialize(fPID);
1106
1107   //------- fcut --------------
1108   fCuts = new AliHFEcuts();
1109   fCuts->CreateStandardCuts();
1110   //fCuts->SetMinNClustersTPC(100);
1111   fCuts->SetMinNClustersTPC(90);
1112   fCuts->SetMinRatioTPCclusters(0.6);
1113   fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
1114   //fCuts->SetMinNClustersITS(3);
1115   fCuts->SetMinNClustersITS(2);
1116   fCuts->SetProductionVertex(0,50,0,50);
1117   fCuts->SetCutITSpixel(AliHFEextraCuts::kAny);
1118   fCuts->SetCheckITSLayerStatus(kFALSE);
1119   fCuts->SetVertexRange(10.);
1120   fCuts->SetTOFPIDStep(kFALSE);
1121   fCuts->SetPtRange(2, 50);
1122   fCuts->SetMaxImpactParam(3.,3.);
1123
1124   //--------Initialize correction Framework and Cuts
1125   fCFM = new AliCFManager;
1126   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
1127   fCFM->SetNStepParticle(kNcutSteps);
1128   for(Int_t istep = 0; istep < kNcutSteps; istep++)
1129     fCFM->SetParticleCutsList(istep, NULL);
1130
1131   if(!fCuts){
1132     AliWarning("Cuts not available. Default cuts will be used");
1133     fCuts = new AliHFEcuts;
1134     fCuts->CreateStandardCuts();
1135   }
1136   fCuts->Initialize(fCFM);
1137
1138   //---------Output Tlist
1139   fOutputList = new TList();
1140   fOutputList->SetOwner();
1141   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
1142
1143   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
1144   fOutputList->Add(fNoEvents);
1145
1146   Int_t binsE[5] =    {250, 100,  1000, 200,   10};
1147   Double_t xminE[5] = {1.0,  -1,   0.0,   0, -0.5};
1148   Double_t xmaxE[5] = {3.5,   1, 100.0, 100,  9.5};
1149   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE);
1150   //if(fqahist==1)fOutputList->Add(fEMCAccE);
1151
1152   hEMCAccE = new TH2F("hEMCAccE","Cluster Energy",200,0,100,100,0,20);
1153   fOutputList->Add(hEMCAccE);
1154
1155   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
1156   fOutputList->Add(fTrkpt);
1157
1158   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
1159   fOutputList->Add(fTrackPtBefTrkCuts);
1160
1161   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
1162   fOutputList->Add(fTrackPtAftTrkCuts);
1163
1164   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
1165   fOutputList->Add(fTPCnsigma);
1166
1167   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
1168   fOutputList->Add(fTrkEovPBef);
1169
1170   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
1171   fOutputList->Add(fTrkEovPAft);
1172
1173   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10);
1174   fOutputList->Add(fdEdxBef);
1175
1176   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10);
1177   fOutputList->Add(fdEdxAft);
1178
1179   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50);
1180   fOutputList->Add(fIncpT);
1181
1182   fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1183   fOutputList->Add(fIncpTM20);
1184
1185   Int_t nBinspho[9] =  { 10,  30,   600, 60,   50,    4,  40,   8,  30};
1186   Double_t minpho[9] = {  0.,  0., -0.1, 40,  0, -0.5,   0,-1.5,   0};
1187   Double_t maxpho[9] = {100., 30.,  0.5, 100,   1,  3.5,   2, 6.5,  30};
1188
1189   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);
1190   if(fqahist==1)fOutputList->Add(fInvmassLS);
1191
1192   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);
1193   if(fqahist==1)fOutputList->Add(fInvmassULS);
1194
1195   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);
1196   if(fqahist==1)fOutputList->Add(fInvmassLSmc);
1197
1198   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);
1199   if(fqahist==1)fOutputList->Add(fInvmassULSmc);
1200
1201   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 );
1202   fInvmassLSreco->Sumw2();
1203   fOutputList->Add(fInvmassLSreco);
1204
1205   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 );
1206   fInvmassULSreco->Sumw2();
1207   fOutputList->Add(fInvmassULSreco);
1208
1209   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 );
1210   fInvmassLSmc0->Sumw2();
1211   fOutputList->Add(fInvmassLSmc0);
1212
1213   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 );
1214   fInvmassLSmc1->Sumw2();
1215   fOutputList->Add(fInvmassLSmc1);
1216
1217   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 );
1218   fInvmassLSmc2->Sumw2();
1219   fOutputList->Add(fInvmassLSmc2);
1220
1221   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 );
1222   fInvmassLSmc3->Sumw2();
1223   fOutputList->Add(fInvmassLSmc3);
1224
1225   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 );
1226   fInvmassULSmc0->Sumw2();
1227   fOutputList->Add(fInvmassULSmc0);
1228
1229   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 );
1230   fInvmassULSmc1->Sumw2();
1231   fOutputList->Add(fInvmassULSmc1);
1232
1233   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 );
1234   fInvmassULSmc2->Sumw2();
1235   fOutputList->Add(fInvmassULSmc2);
1236
1237   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 );
1238   fInvmassULSmc3->Sumw2();
1239   fOutputList->Add(fInvmassULSmc3);
1240
1241   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1242   fOutputList->Add(fOpeningAngleLS);
1243
1244   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1245   fOutputList->Add(fOpeningAngleULS);
1246
1247   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
1248   fOutputList->Add(fPhotoElecPt);
1249
1250   fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50);
1251   fOutputList->Add(fPhoElecPt);
1252
1253   fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50);
1254   fOutputList->Add(fPhoElecPtM20);
1255
1256   fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50);
1257   fOutputList->Add(fSameElecPt);
1258
1259   fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50);
1260   fOutputList->Add(fSameElecPtM20);
1261
1262   fSemiElecPtM20 = new TH2F("fSemiElecPtM20", "Semi-inclusive electron pt with M20",200,0,100,100,0,50);
1263   fOutputList->Add(fSemiElecPtM20);
1264
1265   fCent = new TH1F("fCent","Centrality",200,0,100) ;
1266   fOutputList->Add(fCent);
1267
1268   // Make common binning
1269   const Double_t kMinP = 0.;
1270   const Double_t kMaxP = 20.;
1271
1272   //+++ 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
1273   // 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)
1274   Int_t nBins[16] =  {  100,     7,  60,    20,    90,  100,   25,     60,    60,  100,    40,  10,  250,  100, 100,    8};
1275   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};
1276   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};
1277   //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);
1278   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);
1279   if(fqahist==1)fOutputList->Add(fEleInfo);
1280
1281   // Make common binning
1282   Int_t nBinsEop[3] =  { 10, 50, 100};
1283   Double_t minEop[3] = {  0,  0,  -5};
1284   Double_t maxEop[3] = {100, 50,   5};
1285   fElenSigma= new THnSparseD("fElenSigma", "Electron nSigma; cent; pT [GeV/c]; nSigma", 3, nBinsEop, minEop, maxEop);
1286   fOutputList->Add(fElenSigma);
1287
1288
1289   //<---  Trigger info
1290   /*
1291   fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100);
1292   fOutputList->Add(fClsEBftTrigCut);
1293
1294   fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100);
1295   fOutputList->Add(fClsEAftTrigCut);
1296
1297   fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100);
1298   fOutputList->Add(fClsEAftTrigCut1);
1299
1300   fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100);
1301   fOutputList->Add(fClsEAftTrigCut2);
1302
1303   fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100);
1304   fOutputList->Add(fClsEAftTrigCut3);
1305
1306   fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100);
1307   fOutputList->Add(fClsEAftTrigCut4);
1308
1309   fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1310   fOutputList->Add(fClsETime);
1311
1312   fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009);
1313   fOutputList->Add(fClsETime1);
1314
1315   fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25);
1316   fOutputList->Add(fTrigTimes);
1317
1318   fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000);
1319   fOutputList->Add(fCellCheck);
1320   */
1321   //<---------- MC
1322
1323   fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1324   fOutputList->Add(fInputHFEMC);
1325
1326   fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50);
1327   fOutputList->Add(fInputAlle);
1328
1329   fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50);
1330   fOutputList->Add(fIncpTMChfe);
1331
1332   fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50);
1333   fOutputList->Add(fIncpTMChfeAll);
1334
1335   fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50);
1336   fOutputList->Add(fIncpTMCM20hfe);
1337
1338   fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50);
1339   fOutputList->Add(fIncpTMCM20hfeAll);
1340
1341   fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50);
1342   fOutputList->Add(fIncpTMCM20hfeCheck);
1343
1344   Int_t nBinspho2[5] =  { 200, 100,    7,    3, 700};
1345   Double_t minpho2[5] = {  0.,  0., -2.5, -0.5, 0.};
1346   Double_t maxpho2[5] = {100., 50.,  4.5,  2.5, 70.};
1347
1348   fInputHFEMC_weight = new THnSparseD("fInputHFEMC_weight", "MC HFE electron pt",5,nBinspho2,minpho2,maxpho2);
1349   fOutputList->Add(fInputHFEMC_weight);
1350
1351   fIncpTMCM20hfeCheck_weight = new THnSparseD("fIncpTMCM20hfeCheck_weight", "HFE electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1352   fOutputList->Add(fIncpTMCM20hfeCheck_weight);
1353
1354   fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2);
1355   fOutputList->Add(fIncpTMCpho);
1356
1357   fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1358   fOutputList->Add(fIncpTMCM20pho);
1359
1360   fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1361   fOutputList->Add(fPhoElecPtMC);
1362
1363   fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1364   fOutputList->Add(fPhoElecPtMCM20);
1365
1366   fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2);
1367   fOutputList->Add(fSameElecPtMC);
1368
1369   fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2);
1370   fOutputList->Add(fSameElecPtMCM20);
1371
1372   fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1373   fIncpTMCM20pho_pi0e->Sumw2();
1374   fOutputList->Add(fIncpTMCM20pho_pi0e);
1375
1376   fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1377   fPhoElecPtMCM20_pi0e->Sumw2();
1378   fOutputList->Add(fPhoElecPtMCM20_pi0e);
1379
1380   fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1381   fSameElecPtMCM20_pi0e->Sumw2();
1382   fOutputList->Add(fSameElecPtMCM20_pi0e);
1383  //
1384   fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1385   fIncpTMCM20pho_eta->Sumw2();
1386   fOutputList->Add(fIncpTMCM20pho_eta);
1387
1388   fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2);
1389   fPhoElecPtMCM20_eta->Sumw2();
1390   fOutputList->Add(fPhoElecPtMCM20_eta);
1391
1392   fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1393   fSameElecPtMCM20_eta->Sumw2();
1394   fOutputList->Add(fSameElecPtMCM20_eta);
1395   // ------------
1396   fIncpTMCpho_pi0e_TPC = new THnSparseD("fIncpTMCpho_pi0e_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1397   fIncpTMCpho_pi0e_TPC->Sumw2();
1398   fOutputList->Add(fIncpTMCpho_pi0e_TPC);
1399
1400   fPhoElecPtMC_pi0e_TPC = new THnSparseD("fPhoElecPtMC_pi0e_TPC", "MC Pho-inclusive electron pt with  pi0->e",5,nBinspho2,minpho2,maxpho2);
1401   fPhoElecPtMC_pi0e_TPC->Sumw2();
1402   fOutputList->Add(fPhoElecPtMC_pi0e_TPC);
1403
1404   fSameElecPtMC_pi0e_TPC = new THnSparseD("fSameElecPtMC_pi0e_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1405   fSameElecPtMC_pi0e_TPC->Sumw2();
1406   fOutputList->Add(fSameElecPtMC_pi0e_TPC);
1407
1408   fIncpTMCpho_eta_TPC = new THnSparseD("fIncpTMCpho_eta_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2);
1409   fIncpTMCpho_eta_TPC->Sumw2();
1410   fOutputList->Add(fIncpTMCpho_eta_TPC);
1411
1412   fPhoElecPtMC_eta_TPC = new THnSparseD("fPhoElecPtMC_eta_TPC", "MC Pho-inclusive electron pt with  pi0->e",5,nBinspho2,minpho2,maxpho2);
1413   fPhoElecPtMC_eta_TPC->Sumw2();
1414   fOutputList->Add(fPhoElecPtMC_eta_TPC);
1415
1416   fSameElecPtMC_eta_TPC = new THnSparseD("fSameElecPtMC_eta_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2);
1417   fSameElecPtMC_eta_TPC->Sumw2();
1418   fOutputList->Add(fSameElecPtMC_eta_TPC);
1419
1420
1421   //-------------
1422
1423   CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200);
1424   fOutputList->Add(CheckNclust);
1425
1426   CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5);
1427   fOutputList->Add(CheckNits);
1428
1429   CheckDCA = new TH2D("CheckDCA","DCA check",200,-5,5,200,-5,5);
1430   fOutputList->Add(CheckDCA);
1431   /*
1432   Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5);
1433   fOutputList->Add(Hpi0pTcheck);
1434
1435   HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5);
1436   fOutputList->Add(HETApTcheck);
1437   */
1438
1439   Int_t nBinspho3[3] =  { 200, 100, 3};
1440   Double_t minpho3[3] = {  0.,  0., -0.5};
1441   Double_t maxpho3[3] = {100., 50., 2.5};
1442
1443   Hpi0pTcheck = new THnSparseD("Hpi0pTcheck","Pi0 pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1444   fOutputList->Add(Hpi0pTcheck);
1445
1446   HETApTcheck = new THnSparseD("HETApTcheck","Eta pT from Hijing",3,nBinspho3,minpho3,maxpho3);
1447   fOutputList->Add(HETApTcheck);
1448   //--
1449   HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5);
1450   fOutputList->Add(HphopTcheck);
1451   //
1452   HDpTcheck = new TH2D("HDpTcheck","D pT from Hijing",100,0,50,3,-0.5,2.5);
1453   fOutputList->Add(HDpTcheck);
1454
1455   HBpTcheck = new TH2D("HBpTcheck","B pT from Hijing",100,0,50,3,-0.5,2.5);
1456   fOutputList->Add(HBpTcheck);
1457   //
1458
1459   fpTCheck = new TH1D("fpTCheck","pT check",500,0,50);
1460   fOutputList->Add(fpTCheck);
1461
1462   fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40);
1463   fOutputList->Add(fMomDtoE);
1464
1465   fLabelCheck = new TH2D("fLabelCheck","MC label",50,0,50,5,-1.5,3.5);
1466   fOutputList->Add(fLabelCheck);
1467
1468   fgeoFake = new TH2D("fgeoFake","Label==0 eta and phi",628,0,6.28,200,-1,1);
1469   fOutputList->Add(fgeoFake);
1470
1471   fFakeTrk0 = new TH2D("fFakeTrk0","fake trakcs",10,0,100,20,0,20);
1472   fOutputList->Add(fFakeTrk0);
1473
1474   fFakeTrk1 = new TH2D("fFakeTrk1","true all e a.f. eID",10,0,100,20,0,20);
1475   fOutputList->Add(fFakeTrk1);
1476
1477   ftimingEle = new TH2D("ftimingEle","electron TOF",100,0,20,100,1e-7,1e-6);
1478   fOutputList->Add(ftimingEle);
1479
1480   // eta correction
1481   // note: parameters 01/31new.TPCnSigmaEtaDep
1482   // 70-90 delta_eta = 0.2
1483
1484   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};
1485   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)
1486   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)
1487   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)
1488   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)
1489   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)
1490   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)
1491   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)
1492
1493   fnSigEtaCorr[0] = new TGraphErrors(12,etaval,corr0); // 0-10
1494   fnSigEtaCorr[1] = new TGraphErrors(12,etaval,corr1); // 10-20
1495   fnSigEtaCorr[2] = new TGraphErrors(12,etaval,corr2); // 20-30
1496   fnSigEtaCorr[3] = new TGraphErrors(12,etaval,corr3); // 30-40
1497   fnSigEtaCorr[4] = new TGraphErrors(12,etaval,corr4); // 40-50
1498   fnSigEtaCorr[5] = new TGraphErrors(12,etaval,corr5); // 50-70
1499   fnSigEtaCorr[6] = new TGraphErrors(12,etaval,corr6); // 70-90
1500
1501   fIncMaxE = new TH2D("fIncMaxE","Inc",10,0,100,10,0,100);
1502   fOutputList->Add(fIncMaxE);
1503
1504   fIncReco = new TH2D("fIncReco","Inc",10,0,100,100,0,500);
1505   fOutputList->Add(fIncReco);
1506
1507   fPhoReco = new TH2D("fPhoReco","Pho",10,0,100,100,0,500);
1508   fOutputList->Add(fPhoReco);
1509
1510   fSamReco = new TH2D("fSamReco","Same",10,0,100,100,0,500);
1511   fOutputList->Add(fSamReco);
1512
1513   fIncRecoMaxE = new TH2D("fIncRecoMaxE","Inc",10,0,100,100,0,500);
1514   fOutputList->Add(fIncRecoMaxE);
1515
1516   fPhoRecoMaxE = new TH2D("fPhoRecoMaxE","Pho",10,0,100,100,0,500);
1517   fOutputList->Add(fPhoRecoMaxE);
1518
1519   fSamRecoMaxE = new TH2D("fSamRecoMaxE","Same",10,0,100,100,0,500);
1520   fOutputList->Add(fSamRecoMaxE);
1521
1522   fPhoVertexReco_HFE = new TH2D("fPhoVertexReco_HFE","photon production Vertex mass selection",40,0,20,250,0,50);
1523   fOutputList->Add(fPhoVertexReco_HFE);
1524
1525   fPhoVertexReco_EMCal = new TH2D("fPhoVertexReco_EMCal","photon production Vertex mass selection",40,0,20,250,0,50);
1526   fOutputList->Add(fPhoVertexReco_EMCal);
1527
1528   fPhoVertexReco_Invmass = new TH2D("fPhoVertexReco_Invmass","photon production Vertex mass selection",40,0,20,250,0,50);
1529   fOutputList->Add(fPhoVertexReco_Invmass);
1530
1531   fPhoVertexReco_step0= new TH2D("fPhoVertexReco_step0","photon production Vertex mass selection",40,0,20,250,0,50);
1532   fOutputList->Add(fPhoVertexReco_step0);
1533
1534   fPhoVertexReco_step1= new TH2D("fPhoVertexReco_step1","photon production Vertex mass selection",40,0,20,250,0,50);
1535   fOutputList->Add(fPhoVertexReco_step1);
1536
1537   fPhoVertexReco_step2= new TH2D("fPhoVertexReco_step2","photon production Vertex mass selection",40,0,20,250,0,50);
1538   fOutputList->Add(fPhoVertexReco_step2);
1539
1540   fPhoVertexReco_step3= new TH2D("fPhoVertexReco_step3","photon production Vertex mass selection",40,0,20,250,0,50);
1541   fOutputList->Add(fPhoVertexReco_step3);
1542
1543   fPhoVertexReco_step4= new TH2D("fPhoVertexReco_step4","photon production Vertex mass selection",40,0,20,250,0,50);
1544   fOutputList->Add(fPhoVertexReco_step4);
1545
1546   fPhoVertexReco_step5= new TH2D("fPhoVertexReco_step5","photon production Vertex mass selection",40,0,20,250,0,50);
1547   fOutputList->Add(fPhoVertexReco_step5);
1548
1549   PostData(1,fOutputList);
1550 }
1551
1552 //________________________________________________________________________
1553 void AliAnalysisTaskHFECal::Terminate(Option_t *)
1554 {
1555   // Info("Terminate");
1556         AliAnalysisTaskSE::Terminate();
1557 }
1558
1559 //________________________________________________________________________
1560 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1561 {
1562   // Check single track cuts for a given cut step
1563   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1564   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1565   return kTRUE;
1566 }
1567 //_________________________________________
1568 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)
1569 {
1570   //Identify non-heavy flavour electrons using Invariant mass method
1571
1572   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1573   fTrackCuts->SetRequireTPCRefit(kTRUE);
1574   fTrackCuts->SetRequireITSRefit(kTRUE);
1575   fTrackCuts->SetEtaRange(-0.9,0.9);
1576   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1577   fTrackCuts->SetMinNClustersTPC(90);
1578
1579   //const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
1580   Double_t bfield = fESD->GetMagneticField();
1581   Double_t emass = 0.51*0.001; // (0.51 MeV)
1582
1583   double ptEle = track->Pt();  //add
1584   if(ibgevent==0 && w > 0.0)
1585      {
1586       fpTCheck->Fill(ptEle,w);
1587      }
1588
1589   Bool_t flagPhotonicElec = kFALSE;
1590   Bool_t flagConvinatElec = kFALSE;
1591
1592   int p1 = 0;
1593   if(mce==3)
1594      {
1595        Int_t label = TMath::Abs(track->GetLabel());
1596        TParticle* particle = stack->Particle(label);
1597        p1 = particle->GetFirstMother();
1598      }
1599
1600   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1601   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1602     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1603     if (!trackAsso) {
1604       printf("ERROR: Could not receive track %d\n", jTracks);
1605       continue;
1606     }
1607     if(itrack==jTracks)continue;
1608     int jbgevent = 0;
1609
1610     int p2 = 0;
1611     if(mce==3)
1612     {
1613       Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1614       TParticle* particle2 = stack->Particle(label2);
1615       Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1616       if(MChijing_ass)jbgevent =1;
1617       if(particle2->GetFirstMother()>-1)
1618          p2 = particle2->GetFirstMother();
1619     }
1620
1621     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1622     Double_t mass=999., width = -999;
1623     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1624
1625     //ptPrim = track->Pt();
1626     ptPrim = ptEle;
1627
1628     dEdxAsso = trackAsso->GetTPCsignal();
1629     ptAsso = trackAsso->Pt();
1630     Int_t chargeAsso = trackAsso->Charge();
1631     Int_t charge = track->Charge();
1632
1633
1634     if(ptAsso <fMimpTassCut) continue;
1635     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1636     //if(dEdxAsso <65 || dEdxAsso>100) continue;
1637     double fTPCnSigmaAss = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
1638     //cout << "fTPCnSigmaAss = " << fTPCnSigmaAss << endl;
1639     //cout << "fTPCnSigmaAss Cut = " << fMimNsigassCut << endl;
1640     if(fTPCnSigmaAss <fMimNsigassCut || fTPCnSigmaAss>5) continue;
1641     //cout << "fTPCnSigmaAss a.f. cut = " << fTPCnSigmaAss << endl;
1642
1643     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1644     if(charge>0) fPDGe1 = -11;
1645     if(chargeAsso>0) fPDGe2 = -11;
1646
1647     //printf("chargeAsso = %d\n",chargeAsso);
1648     //printf("charge = %d\n",charge);
1649     if(charge == chargeAsso) fFlagLS = kTRUE;
1650     if(charge != chargeAsso) fFlagULS = kTRUE;
1651
1652     //printf("fFlagLS = %d\n",fFlagLS);
1653     //printf("fFlagULS = %d\n",fFlagULS);
1654     printf("\n");
1655
1656     AliKFParticle::SetField(bfield);
1657     AliKFParticle ge1(*track, fPDGe1);
1658     AliKFParticle ge2(*trackAsso, fPDGe2);
1659
1660     if(fSetMassNonlinear)
1661       {
1662        ge1.SetNonlinearMassConstraint(emass);
1663        ge2.SetNonlinearMassConstraint(emass);
1664       }
1665
1666     AliKFParticle recg(ge1, ge2);
1667
1668     /*
1669     // vertex
1670     AliKFVertex primV(*pVtx);
1671     primV += recg;
1672     primV -= ge1;
1673     primV -= ge2;
1674     recg.SetProductionVertex(primV);
1675     */
1676
1677     // check chi2
1678     if(recg.GetNDF()<1) continue;
1679     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1680     Double_t chi2cut = 3.0;
1681
1682     // mass.
1683     if(fSetMassConstraint)
1684       {
1685        recg.SetMassConstraint(0,0.0001);
1686        chi2cut = 30.0;
1687       }
1688     recg.GetMass(mass,width);
1689
1690     if(fSetMassWidthCut && width>1e10)continue;
1691
1692         // angle
1693     openingAngle = ge1.GetAngle(ge2);
1694     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1695     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1696
1697     double ishower = 0;
1698     if(shower>0.0 && shower<0.3)ishower = 1;
1699
1700     double phoinfo[9];
1701     phoinfo[0] = cent;
1702     phoinfo[1] = ptPrim;
1703     phoinfo[2] = mass;
1704     phoinfo[3] = nSig;
1705     //phoinfo[3] = dEdxAsso;
1706     phoinfo[4] = openingAngle;
1707     phoinfo[5] = ishower;
1708     phoinfo[6] = ep;
1709     phoinfo[7] = mce;
1710     phoinfo[8] = ptAsso;
1711
1712     if(fFlagLS) fInvmassLS->Fill(phoinfo);
1713     if(fFlagULS) fInvmassULS->Fill(phoinfo);
1714     if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1715     if(fFlagULS && ibgevent==0 && jbgevent==0)
1716        {
1717          fInvmassULSmc->Fill(phoinfo,w);
1718        }
1719     //printf("fInvmassCut %f\n",fInvmassCut);
1720     //printf("openingAngle %f\n",fOpeningAngleCut);
1721
1722     // angle cut
1723     if(openingAngle > fOpeningAngleCut) continue;
1724     // chi2 cut
1725     //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue;
1726     if(chi2recg>chi2cut) continue;
1727
1728     if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass);
1729     if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass);
1730
1731     // for real data
1732     //printf("mce =%f\n",mce);
1733     if(mce<-0.5) // mce==-1. is real
1734        {
1735          //printf("Real data\n");
1736          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1737                flagPhotonicElec = kTRUE;
1738               }
1739          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1740                flagConvinatElec = kTRUE;
1741               }
1742         }
1743     // for MC data
1744     else
1745        {
1746          //printf("MC data\n");
1747
1748          if(w>0.0)
1749            {
1750            //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1751            if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass);
1752            if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass);
1753            if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass);
1754            if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass);
1755            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass);
1756            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass);
1757            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass);
1758            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass);
1759           }
1760
1761          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1762          //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ //<--- only MC train (55,56) v5-03-68-AN , 69 & v5-05-70-AN (till 74)
1763                flagPhotonicElec = kTRUE;
1764               }
1765          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1766                flagConvinatElec = kTRUE;
1767               }
1768         }
1769
1770   }
1771   fFlagPhotonicElec = flagPhotonicElec;
1772   fFlagConvinatElec = flagConvinatElec;
1773
1774 }
1775
1776 //_________________________________________
1777 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)
1778 {
1779   //Identify non-heavy flavour electrons using Invariant mass method
1780
1781   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
1782   fTrackCuts->SetRequireTPCRefit(kTRUE);
1783   fTrackCuts->SetRequireITSRefit(kTRUE);
1784   fTrackCuts->SetEtaRange(-0.9,0.9);
1785   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
1786   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
1787   fTrackCuts->SetMinNClustersTPC(90);
1788
1789   Double_t eMass = TDatabasePDG::Instance()->GetParticle(11)->Mass(); //Electron mass in GeV
1790   Double_t bfield = fESD->GetMagneticField();
1791
1792   double ptEle = track->Pt();  //add
1793   if(ibgevent==0 && w > 0.0)
1794      {
1795       fpTCheck->Fill(ptEle,w);
1796      }
1797
1798   Bool_t flagPhotonicElec = kFALSE;
1799   Bool_t flagConvinatElec = kFALSE;
1800
1801   int p1 = 0;
1802   if(mce==3)
1803      {
1804        Int_t label = TMath::Abs(track->GetLabel());
1805        TParticle* particle = stack->Particle(label);
1806        p1 = particle->GetFirstMother();
1807      }
1808
1809   //for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1810   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
1811     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
1812     if (!trackAsso) {
1813       printf("ERROR: Could not receive track %d\n", jTracks);
1814       continue;
1815     }
1816     if(itrack==jTracks)continue;
1817     int jbgevent = 0;
1818
1819     int p2 = 0;
1820     if(mce==3)
1821     {
1822       Int_t label2 = TMath::Abs(trackAsso->GetLabel());
1823       TParticle* particle2 = stack->Particle(label2);
1824       Bool_t MChijing_ass = fMC->IsFromBGEvent(label2);
1825       if(MChijing_ass)jbgevent =1;
1826       if(particle2->GetFirstMother()>-1)
1827          p2 = particle2->GetFirstMother();
1828     }
1829
1830     Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.;
1831     //Double_t mass=999., width = -999;
1832     Double_t mass=999.;
1833     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1834
1835     //ptPrim = track->Pt();
1836     ptPrim = ptEle;
1837
1838     dEdxAsso = trackAsso->GetTPCsignal();
1839     ptAsso = trackAsso->Pt();
1840     Int_t chargeAsso = trackAsso->Charge();
1841     Int_t charge = track->Charge();
1842
1843
1844     if(ptAsso <0.5) continue;
1845     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
1846     if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1
1847
1848     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1849     if(charge>0) fPDGe1 = -11;
1850     if(chargeAsso>0) fPDGe2 = -11;
1851
1852     //printf("chargeAsso = %d\n",chargeAsso);
1853     //printf("charge = %d\n",charge);
1854     if(charge == chargeAsso) fFlagLS = kTRUE;
1855     if(charge != chargeAsso) fFlagULS = kTRUE;
1856
1857     // mass cal 0
1858     double xt1 = 0.0, xt2 = 0.0;
1859     Double_t p1at[3] = {0,0,0};
1860     Double_t p2at[3] = {0,0,0};
1861     //double dca12 = trackAsso->GetDCA(track,bfield,xt2,xt1);           //DCA track1-track2
1862     double kHasdcaT1 = track->GetPxPyPzAt(xt1,bfield,p1at);             //Track1
1863     double kHasdcaT2 = trackAsso->GetPxPyPzAt(xt2,bfield,p2at);         //Track2
1864    if(!kHasdcaT1 || !kHasdcaT2) AliWarning("It could be a problem in the extrapolation");
1865    //cout << "dca = " << dca12 << endl;
1866    // 3D
1867    TLorentzVector electron1;
1868    TLorentzVector electron2;
1869    TLorentzVector mother;
1870
1871    electron1.SetXYZM(p1at[0], p1at[1], p1at[2], eMass);
1872    electron2.SetXYZM(p2at[0], p2at[1], p2at[2], eMass);
1873    openingAngle = TVector2::Phi_0_2pi(electron1.Angle(electron2.Vect()));
1874
1875    mother      = electron1 + electron2;
1876    //double invmassAtDCA  = mother.M();
1877
1878    // 2D
1879    TLorentzVector electron1_2D;
1880    TLorentzVector electron2_2D;
1881    TLorentzVector mother_2D;
1882
1883    double pT1at = sqrt(pow(p1at[0],2)+pow(p1at[1],2));
1884    double pT2at = sqrt(pow(p2at[0],2)+pow(p2at[1],2));
1885
1886    electron1_2D.SetXYZM(pT1at, 0.0, p1at[2], eMass);
1887    electron2_2D.SetXYZM(pT2at, 0.0, p2at[2], eMass);
1888
1889    mother_2D      = electron1_2D + electron2_2D;
1890    double invmassAtDCA_2D  = mother_2D.M();
1891    mass = invmassAtDCA_2D;
1892
1893         // angle
1894     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
1895     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
1896
1897     double ishower = 0;
1898     if(shower>0.0 && shower<0.3)ishower = 1;
1899
1900     double phoinfo[9];
1901     phoinfo[0] = cent;
1902     phoinfo[1] = ptPrim;
1903     phoinfo[2] = mass;
1904     phoinfo[3] = nSig;
1905     //phoinfo[3] = dEdxAsso;
1906     phoinfo[4] = openingAngle;
1907     phoinfo[5] = ishower;
1908     phoinfo[6] = ep;
1909     phoinfo[7] = mce;
1910     phoinfo[8] = ptAsso;
1911
1912     if(fFlagLS) fInvmassLS->Fill(phoinfo);
1913     if(fFlagULS) fInvmassULS->Fill(phoinfo);
1914     if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w);
1915     if(fFlagULS && ibgevent==0 && jbgevent==0)
1916        {
1917          fInvmassULSmc->Fill(phoinfo,w);
1918        }
1919     //printf("fInvmassCut %f\n",fInvmassCut);
1920     //printf("openingAngle %f\n",fOpeningAngleCut);
1921
1922     // angle cut
1923     if(openingAngle > fOpeningAngleCut) continue;
1924     // chi2 cut
1925     //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue;
1926     //if(chi2recg>chi2cut) continue;
1927
1928     if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass);
1929     if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass);
1930
1931     // for real data
1932     //printf("mce =%f\n",mce);
1933     if(mce<-0.5) // mce==-1. is real
1934        {
1935          //printf("Real data\n");
1936          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
1937                flagPhotonicElec = kTRUE;
1938               }
1939          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec){
1940                flagConvinatElec = kTRUE;
1941               }
1942         }
1943     // for MC data
1944     else
1945        {
1946          //printf("MC data\n");
1947
1948          if(w>0.0)
1949            {
1950            //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl;
1951            if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass);
1952            if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass);
1953            if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass);
1954            if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass);
1955            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass);
1956            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass);
1957            if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass);
1958            if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass);
1959           }
1960
1961          if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (ibgevent==jbgevent)){
1962          //if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec && (p1==p2)){ <--- only MC train (55,56) v5-03-68-AN & 69 for check
1963                flagPhotonicElec = kTRUE;
1964               }
1965          if(mass<fInvmassCut && fFlagLS && !flagConvinatElec && (ibgevent==jbgevent)){
1966                flagConvinatElec = kTRUE;
1967               }
1968         }
1969
1970   }
1971   fFlagPhotonicElec = flagPhotonicElec;
1972   fFlagConvinatElec = flagConvinatElec;
1973
1974 }
1975
1976 //-------------------------------------------
1977
1978 void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid)
1979 {
1980  //int label = 99999;
1981  //int pid = 99999;
1982
1983  if(part->GetFirstMother()>-1)
1984    {
1985     label = part->GetFirstMother();
1986     pid = stack->Particle(label)->GetPdgCode();
1987    }
1988    //cout << "Find Mother : label = " << label << " ; pid" << pid << endl;
1989 }
1990
1991 double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT)
1992 {
1993         double weight = 1.0;
1994
1995         if(mcPi0pT>0.0 && mcPi0pT<5.0)
1996         {
1997                 weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT));
1998         }
1999         else
2000         {
2001                 weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65));
2002         }
2003   return weight;
2004 }
2005
2006 double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT)
2007 {
2008   double weight = 1.0;
2009
2010   weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65);
2011   return weight;
2012 }
2013
2014
2015 //_________________________________________
2016 void AliAnalysisTaskHFECal::FindTriggerClusters()
2017 {
2018   //cout << "finding trigger patch" << endl;
2019   // constants
2020   const int nModuleCols = 2;
2021   const int nModuleRows = 5;
2022   const int nColsFeeModule = 48;
2023   const int nRowsFeeModule = 24;
2024   const int nColsFaltroModule = 24;
2025   const int nRowsFaltroModule = 12;
2026   //const int faltroWidthMax = 20;
2027
2028   // part 1, trigger extraction -------------------------------------
2029   Int_t globCol, globRow;
2030   //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0;
2031   Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0;
2032
2033   //Int_t trigtimes[faltroWidthMax];
2034   Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
2035   Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows];
2036   //Double_t fTrigCutLow = 6;
2037   //Double_t fTrigCutHigh = 10;
2038   Double_t fTimeCutLow =  469e-09;
2039   Double_t fTimeCutHigh = 715e-09;
2040
2041   AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" );
2042
2043   // erase trigger maps
2044   for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ )
2045   {
2046     for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ )
2047     {
2048       ftriggersCut[i][j] = 0;
2049       ftriggers[i][j] = 0;
2050       ftriggersTime[i][j] = 0;
2051     }
2052   }
2053
2054   Int_t iglobCol=0, iglobRow=0;
2055   // go through triggers
2056   if( fCaloTrigger->GetEntries() > 0 )
2057   {
2058     // needs reset
2059     fCaloTrigger->Reset();
2060     while( fCaloTrigger->Next() )
2061     {
2062       fCaloTrigger->GetPosition( globCol, globRow );
2063       fCaloTrigger->GetNL0Times( ntimes );
2064       /*
2065       // no L0s
2066       if( ntimes < 1 )   continue;
2067       // get precise timings
2068       fCaloTrigger->GetL0Times( trigtimes );
2069       trigInCut = 0;
2070       for(Int_t i = 0; i < ntimes; i++ )
2071       {
2072         // save the first trigger time in channel
2073         if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] )
2074           triggersTime[globCol][globRow] = trigtimes[i];
2075         //printf("trigger times: %d\n",trigtimes[i]);
2076         // check if in cut
2077         if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh )
2078           trigInCut = 1;
2079
2080         fTrigTimes->Fill(trigtimes[i]);
2081       }
2082       */
2083
2084       //L1 analysis from AliAnalysisTaskEMCALTriggerQA
2085       Int_t bit = 0;
2086       fCaloTrigger->GetTriggerBits(bit);
2087       //cout << "bit = " << bit << endl;
2088
2089       Int_t ts = 0;
2090       fCaloTrigger->GetL1TimeSum(ts);
2091       //cout << "ts = " << ts << endl;
2092       if (ts > 0)ftriggers[globCol][globRow] = 1;
2093       // number of triggered channels in event
2094       nTrigChannel++;
2095       // ... inside cut
2096       if(ts>0 && (bit >> 6 & 0x1))
2097       {
2098         iglobCol = globCol;
2099         iglobRow = globRow;
2100         nTrigChannelCut++;
2101         //cout << "ts cut = " << ts << endl;
2102         //cout << "globCol = " << globCol << endl;
2103         //cout << "globRow = " << globRow << endl;
2104         ftriggersCut[globCol][globRow] = 1;
2105       }
2106
2107     } // calo trigger entries
2108   } // has calo trigger entries
2109
2110   // part 2 go through the clusters here -----------------------------------
2111   //cout << " part 2 go through the clusters here ----------------------------------- " << endl;
2112   Int_t nCluster=0, nCell=0, iCell=0, gCell=0;
2113   Short_t cellAddr, nSACell;
2114   Int_t mclabel;
2115   //Int_t nSACell, iSACell, mclabel;
2116   Int_t iSACell;
2117   Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0;
2118   Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi;
2119
2120   TRefArray *fCaloClusters = new TRefArray();
2121   fESD->GetEMCALClusters( fCaloClusters );
2122   nCluster = fCaloClusters->GetEntries();
2123
2124
2125   // save all cells times if there are clusters
2126   if( nCluster > 0 ){
2127     // erase time map
2128     for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){
2129       for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){
2130         cellTime[i][j] = 0.;
2131         cellEnergy[i][j] = 0.;
2132       }
2133     }
2134
2135     // get cells
2136     AliESDCaloCells *fCaloCells = fESD->GetEMCALCells();
2137     //AliVCaloCells fCaloCells = fESD->GetEMCALCells();
2138     nSACell = fCaloCells->GetNumberOfCells();
2139     for(iSACell = 0; iSACell < nSACell; iSACell++ ){
2140       // get the cell info *fCal
2141       fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac);
2142
2143       // get cell position
2144       fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta );
2145       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2146
2147       // convert co global phi eta
2148       gphi = iphi + nRowsFeeModule*(nSupMod/2);
2149       geta = ieta + nColsFeeModule*(nSupMod%2);
2150
2151       // save cell time and energy
2152       cellTime[geta][gphi] = cellTimeT;
2153       cellEnergy[geta][gphi] = cellAmp;
2154
2155     }
2156   }
2157
2158   Int_t nClusterTrig, nClusterTrigCut;
2159   UShort_t *cellAddrs;
2160   Double_t clsE=-999, clsEta=-999, clsPhi=-999;
2161   Float_t clsPos[3] = {0.,0.,0.};
2162
2163   for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++)
2164   {
2165     AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl);
2166     if(!cluster || !cluster->IsEMCAL()) continue;
2167
2168     // get cluster cells
2169     nCell = cluster->GetNCells();
2170
2171     // get cluster energy
2172     clsE = cluster->E();
2173
2174     // get cluster position
2175     cluster->GetPosition(clsPos);
2176     TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
2177     clsEta = clsPosVec.Eta();
2178     clsPhi = clsPosVec.Phi();
2179
2180     // get the cell addresses
2181     cellAddrs = cluster->GetCellsAbsId();
2182
2183     // check if the cluster contains cell, that was marked as triggered
2184     nClusterTrig = 0;
2185     nClusterTrigCut = 0;
2186
2187     // loop the cells to check, if cluser in acceptance
2188     // any cluster with a cell outside acceptance is not considered
2189     for( iCell = 0; iCell < nCell; iCell++ )
2190     {
2191      // check hot cell
2192      //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]);
2193
2194       // get cell position
2195       fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta );
2196       fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta);
2197
2198       // convert co global phi eta
2199       gphi = iphi + nRowsFeeModule*(nSupMod/2);
2200       geta = ieta + nColsFeeModule*(nSupMod%2);
2201
2202       if( cellTime[geta][gphi] > 0. ){
2203         clusterTime += cellTime[geta][gphi];
2204         gCell++;
2205       }
2206
2207       // get corresponding FALTRO
2208       fphi = gphi / 2;
2209       feta = geta / 2;
2210
2211         //cout << "fphi = " << fphi << endl;
2212         //cout << "feta = " << feta << endl;
2213
2214       // try to match with a triggered
2215       if( ftriggers[feta][fphi]==1)
2216       {  nClusterTrig++;
2217       }
2218       if( ftriggersCut[feta][fphi]==1)
2219       { nClusterTrigCut++;
2220       }
2221
2222       //cout << "nClusterTrigCut : " << nClusterTrigCut << endl;
2223
2224     } // cells
2225
2226
2227     if( gCell > 0 )
2228       clusterTime = clusterTime / (Double_t)gCell;
2229     // fix the reconstriction code time 100ns jumps
2230     if( fESD->GetBunchCrossNumber() % 4 < 2 )
2231       clusterTime -= 0.0000001;
2232
2233     //fClsETime->Fill(clsE,clusterTime);
2234     //fClsEBftTrigCut->Fill(clsE);
2235
2236     if(nClusterTrig>0){
2237       //fClsETime1->Fill(clsE,clusterTime);
2238     }
2239
2240     if(nClusterTrig>0){
2241       cluster->SetChi2(1);
2242       //fClsEAftTrigCut1->Fill(clsE);
2243     }
2244
2245     if(nClusterTrigCut>0){
2246       cluster->SetChi2(2);
2247       //fClsEAftTrigCut2->Fill(clsE);
2248     }
2249
2250     if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3)))
2251     {
2252       cluster->SetChi2(3);
2253       //fClsEAftTrigCut3->Fill(clsE);
2254     }
2255
2256     if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh ))
2257     {
2258       // cluster->SetChi2(4);
2259       //fClsEAftTrigCut4->Fill(clsE);
2260     }
2261     if(nClusterTrigCut<1)
2262     {
2263       cluster->SetChi2(0);
2264
2265       //fClsEAftTrigCut->Fill(clsE);
2266     }
2267
2268   } // clusters
2269 }
2270
2271 // <-------- only MC correction
2272 double AliAnalysisTaskHFECal::MCEopMeanCorrection(double pTmc, float central)
2273 {
2274   TF1 *fcorr0 = new TF1("fcorr0","[0]*tanh([1]+[2]*x)");
2275   TF1 *fcorr1 = new TF1("fcorr1","[0]*tanh([1]+[2]*x)");
2276
2277  double shift = 0.0;
2278
2279  if(central>0 && central<=10)
2280    {
2281     fcorr0->SetParameters(1.045,1.288,3.18e-01); //
2282     fcorr1->SetParameters(9.91e-01,3.466,2.344);
2283    }
2284  else if(central>10 && central<=20)
2285    {
2286     fcorr0->SetParameters(1.029,8.254e-01,4.07e-01);
2287     fcorr1->SetParameters(0.975,2.276,1.501e-01);
2288    }
2289  else if(central>20 && central<=30)
2290    {
2291     fcorr0->SetParameters(1.01,8.795e-01,3.904e-01);
2292     fcorr1->SetParameters(9.675e-01,1.654,2.583e-01);
2293    }
2294  else if(central>30 && central<=40)
2295    {
2296     fcorr0->SetParameters(1.00,1.466,2.305e-1);
2297     fcorr1->SetParameters(9.637e-01,1.473,2.754e-01);
2298    }
2299  else if(central>40 && central<=50)
2300    {
2301     fcorr0->SetParameters(1.00,1.422,1.518e-01);
2302     fcorr1->SetParameters(9.59e-01,1.421,2.931e-01);
2303    }
2304
2305  else if(central>50 && central<=70)
2306    {
2307     fcorr0->SetParameters(0.989,2.495,2.167);
2308     fcorr1->SetParameters(0.961,1.734,1.438e-01);
2309    }
2310  else if(central>70 && central<=100)
2311    {
2312     fcorr0->SetParameters(0.981,-3.373,3.93327);
2313     fcorr1->SetParameters(9.574e-01,1.698,1.58e-01);
2314    }
2315
2316
2317  shift = fcorr0->Eval(pTmc)-fcorr1->Eval(pTmc);
2318
2319  fcorr0->Delete();
2320  fcorr1->Delete();
2321
2322  return shift;
2323 }
2324
2325 // <-------- only Data correction
2326 double AliAnalysisTaskHFECal::NsigmaCorrection(double tmpeta, float central)
2327 {
2328  int icent = 0;
2329
2330  if(central>=0 && central<10)
2331    {
2332     icent = 0;
2333    }
2334  else if(central>=10 && central<20)
2335   {
2336    icent = 1;
2337   }
2338  else if(central>=20 && central<30)
2339   {
2340    icent = 2;
2341   }
2342  else if(central>=30 && central<40)
2343   {
2344    icent = 3;
2345   }
2346  else if(central>=40 && central<50)
2347   {
2348    icent = 4;
2349   }
2350  else if(central>=50 && central<70)
2351   {
2352    icent = 5;
2353   }
2354  else
2355   {
2356    icent = 6;
2357   }
2358
2359  double shift = fnSigEtaCorr[icent]->Eval(tmpeta);
2360
2361  //cout << "eta correction"<< endl;
2362  //cout << "cent = "<< central<< endl;
2363  //cout << "icent = "<< icent << endl;
2364  //cout << "shift = "<< shift << endl;
2365
2366  return shift;
2367
2368 }
2369
2370
2371 double AliAnalysisTaskHFECal::SumpT(Int_t itrack, AliESDtrack* track)
2372 {
2373
2374   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
2375   fTrackCuts->SetRequireTPCRefit(kTRUE);
2376   fTrackCuts->SetRequireITSRefit(kTRUE);
2377   fTrackCuts->SetEtaRange(-0.9,0.9);
2378   //fTrackCuts->SetRequireSigmaToVertex(kTRUE);
2379   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
2380   fTrackCuts->SetMinNClustersTPC(90);
2381
2382   double pTrecp = track->Pt();
2383   double phiorg = track->Phi();
2384   double etaorg = track->Eta();
2385
2386   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
2387     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
2388     if (!trackAsso) {
2389       printf("ERROR: Could not receive track %d\n", jTracks);
2390       continue;
2391     }
2392     if(itrack==jTracks)continue;
2393     double pTAss = trackAsso->Pt();
2394     double etaAss = trackAsso->Eta();
2395     double phiAss = trackAsso->Phi();
2396
2397     double delphi = phiorg - phiAss;
2398     double deleta = etaorg - etaAss;
2399
2400     double R = sqrt(pow(deleta,2)+pow(delphi,2));
2401     if(pTAss<0.5)continue;
2402     if(R<0.4)pTrecp+=pTAss;
2403
2404     }
2405
2406    return pTrecp;
2407 }
2408