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