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