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