PWGHF/hfe/AliAnalysisTaskElecV2.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 v2 with EMCal triggered events
  17 // Author: Denise Godoy
  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 "AliAnalysisTaskElecV2.h"
  39 #include "TGeoGlobalMagField.h"
  40 #include "AliLog.h"
  41 #include "AliAnalysisTaskSE.h"
  42 #include "TRefArray.h"
  43 #include "TVector.h"
  44 #include "AliESDInputHandler.h"
  45 #include "AliESDpid.h"
  46 #include "AliESDtrackCuts.h"
  47 #include "AliPhysicsSelection.h"
  48 #include "AliESDCaloCluster.h"
  49 #include "AliAODCaloCluster.h"
  50 #include "AliEMCALRecoUtils.h"
  51 #include "AliEMCALGeometry.h"
  52 #include "AliGeomManager.h"
  53 #include "stdio.h"
  54 #include "TGeoManager.h"
  55 #include "iostream"
  56 #include "fstream"
  57
  58 #include "AliEMCALTrack.h"
  59 #include "AliMagF.h"
  60
  61 #include "AliKFParticle.h"
  62 #include "AliKFVertex.h"
  63
  64 #include "AliMCEventHandler.h"
  65 #include "AliMCEvent.h"
  66 #include "AliMCParticle.h"
  67 #include "AliStack.h"
  68
  69 #include "AliPID.h"
  70 #include "AliPIDResponse.h"
  71 #include "AliHFEcontainer.h"
  72 #include "AliHFEcuts.h"
  73 #include "AliHFEpid.h"
  74 #include "AliHFEpidBase.h"
  75 #include "AliHFEpidQAmanager.h"
  76 #include "AliHFEtools.h"
  77 #include "AliCFContainer.h"
  78 #include "AliCFManager.h"
  79
  80 #include "AliEventplane.h"
  81 #include "AliCentrality.h"
  82
  83 ClassImp(AliAnalysisTaskElecV2)
  84 //________________________________________________________________________
  85 AliAnalysisTaskElecV2::AliAnalysisTaskElecV2(const char *name)
  86   : AliAnalysisTaskSE(name)
  87   ,fESD(0)
  88   ,fMC(0)
  89   ,fOutputList(0)
  90   ,fTrackCuts(0)
  91   ,fCuts(0)
  92   ,fIdentifiedAsOutInz(kFALSE)
  93   ,fPassTheEventCut(kFALSE)
  94   ,fRejectKinkMother(kFALSE)
  95   ,fIsMC(kFALSE)
  96   ,fVz(0.0)
  97   ,fCFM(0)
  98   ,fPID(0)
  99   ,fPIDqa(0)
 100   ,fOpeningAngleCut(0.1)
 101   ,fInvmassCut(0.01)
 102   ,fNoEvents(0)
 103   ,fTrkpt(0)
 104   ,fTrkEovPBef(0)
 105   ,fTrkEovPAft(0)
 106   ,fdEdxBef(0)
 107   ,fdEdxAft(0)
 108   ,fInvmassLS(0)
 109   ,fInvmassULS(0)
 110   ,fOpeningAngleLS(0)
 111   ,fOpeningAngleULS(0)
 112   ,fPhotoElecPt(0)
 113   ,fSemiInclElecPt(0)
 114   ,fMCphotoElecPt(0)
 115   ,fTrackPtBefTrkCuts(0)
 116   ,fTrackPtAftTrkCuts(0)
 117   ,fTPCnsigma(0)
 118   ,fCent(0)
 119   ,fevPlaneV0A(0)
 120   ,fevPlaneV0C(0)
 121   ,fevPlaneTPC(0)
 122   ,fTPCsubEPres(0)
 123   ,fEPres(0)
 124   ,fCorr(0)
 125   ,feTPCV2(0)
 126   ,feV2(0)
 127   ,fphoteV2(0)
 128   ,fChargPartV2(0)
 129   ,fGammaWeight(0)
 130   ,fPi0Weight(0)
 131   ,fEtaWeight(0)
 132 {
 133   //Named constructor
 134
 135   fPID = new AliHFEpid("hfePid");
 136   fTrackCuts = new AliESDtrackCuts();
 137
 138   // Define input and output slots here
 139   // Input slot #0 works with a TChain
 140   DefineInput(0, TChain::Class());
 141   // Output slot #0 id reserved by the base class for AOD
 142   // Output slot #1 writes into a TH1 container
 143   // DefineOutput(1, TH1I::Class());
 144   DefineOutput(1, TList::Class());
 145   //  DefineOutput(3, TTree::Class());
 146 }
 147
 148 //________________________________________________________________________
 149 AliAnalysisTaskElecV2::AliAnalysisTaskElecV2()
 150   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
 151   ,fESD(0)
 152   ,fMC(0)
 153   ,fOutputList(0)
 154   ,fTrackCuts(0)
 155   ,fCuts(0)
 156   ,fIdentifiedAsOutInz(kFALSE)
 157   ,fPassTheEventCut(kFALSE)
 158   ,fRejectKinkMother(kFALSE)
 159   ,fIsMC(kFALSE)
 160   ,fVz(0.0)
 161   ,fCFM(0)
 162   ,fPID(0)
 163   ,fPIDqa(0)
 164   ,fOpeningAngleCut(0.1)
 165   ,fInvmassCut(0.01)
 166   ,fNoEvents(0)
 167   ,fTrkpt(0)
 168   ,fTrkEovPBef(0)
 169   ,fTrkEovPAft(0)
 170   ,fdEdxBef(0)
 171   ,fdEdxAft(0)
 172   ,fInvmassLS(0)
 173   ,fInvmassULS(0)
 174   ,fOpeningAngleLS(0)
 175   ,fOpeningAngleULS(0)
 176   ,fPhotoElecPt(0)
 177   ,fSemiInclElecPt(0)
 178   ,fMCphotoElecPt(0)
 179   ,fTrackPtBefTrkCuts(0)
 180   ,fTrackPtAftTrkCuts(0)
 181   ,fTPCnsigma(0)
 182   ,fCent(0)
 183   ,fevPlaneV0A(0)
 184   ,fevPlaneV0C(0)
 185   ,fevPlaneTPC(0)
 186   ,fTPCsubEPres(0)
 187   ,fEPres(0)
 188   ,fCorr(0)
 189   ,feTPCV2(0)
 190   ,feV2(0)
 191   ,fphoteV2(0)
 192   ,fChargPartV2(0)
 193   ,fGammaWeight(0)
 194   ,fPi0Weight(0)
 195   ,fEtaWeight(0)
 196 {
 197         //Default constructor
 198         fPID = new AliHFEpid("hfePid");
 199
 200         fTrackCuts = new AliESDtrackCuts();
 201
 202         // Constructor
 203         // Define input and output slots here
 204         // Input slot #0 works with a TChain
 205         DefineInput(0, TChain::Class());
 206         // Output slot #0 id reserved by the base class for AOD
 207         // Output slot #1 writes into a TH1 container
 208         // DefineOutput(1, TH1I::Class());
 209         DefineOutput(1, TList::Class());
 210         //DefineOutput(3, TTree::Class());
 211 }
 212 //_________________________________________
 213
 214 AliAnalysisTaskElecV2::~AliAnalysisTaskElecV2()
 215 {
 216   //Destructor
 217
 218   delete fOutputList;
 219   delete fPID;
 220   delete fCFM;
 221   delete fPIDqa;
 222   delete fTrackCuts;
 223 }
 224 //_________________________________________
 225
 226 void AliAnalysisTaskElecV2::UserExec(Option_t*)
 227 {
 228   //Main loop
 229   //Called for each event
 230
 231   // create pointer to event
 232   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 233   if (!fESD) {
 234     printf("ERROR: fESD not available\n");
 235     return;
 236   }
 237
 238   if(!fCuts){
 239     AliError("HFE cuts not available");
 240     return;
 241   }
 242
 243   if(!fPID->IsInitialized()){
 244     // Initialize PID with the given run number
 245     AliWarning("PID not initialised, get from Run no");
 246     fPID->InitializePID(fESD->GetRunNumber());
 247   }
 248
 249   if(fIsMC)fMC = MCEvent();
 250   AliStack* stack = NULL;
 251   if(fIsMC && fMC) stack = fMC->Stack();
 252
 253   Int_t fNOtrks =  fESD->GetNumberOfTracks();
 254   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 255
 256   Double_t pVtxZ = -999;
 257   pVtxZ = pVtx->GetZ();
 258
 259   if(TMath::Abs(pVtxZ)>10) return;
 260   fNoEvents->Fill(0);
 261
 262   if(fNOtrks<2) return;
 263
 264   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 265   if(!pidResponse){
 266     AliDebug(1, "Using default PID Response");
 267     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 268   }
 269
 270   fPID->SetPIDResponse(pidResponse);
 271
 272   fCFM->SetRecEventInfo(fESD);
 273
 274   Float_t cent = -1.;
 275   AliCentrality *centrality = fESD->GetCentrality();
 276   cent = centrality->GetCentralityPercentile("V0M");
 277   fCent->Fill(cent);
 278
 279   //Event planes
 280
 281   Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
 282   if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
 283   fevPlaneV0A->Fill(evPlaneV0A,cent);
 284
 285   Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
 286   if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
 287   fevPlaneV0C->Fill(evPlaneV0C,cent);
 288
 289   AliEventplane* esdTPCep = fESD->GetEventplane();
 290   TVector2 *standardQ = 0x0;
 291   Double_t qx = -999., qy = -999.;
 292   standardQ = esdTPCep->GetQVector();
 293   if(!standardQ)return;
 294
 295   qx = standardQ->X();
 296   qy = standardQ->Y();
 297
 298   TVector2 qVectorfortrack;
 299   qVectorfortrack.Set(qx,qy);
 300   Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
 301   fevPlaneTPC->Fill(evPlaneTPC,cent);
 302
 303   TVector2 *qsub1a = esdTPCep->GetQsub1();
 304   TVector2 *qsub2a = esdTPCep->GetQsub2();
 305   Double_t evPlaneResTPC = -999.;
 306   if(qsub1a && qsub2a)
 307   {
 308           evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
 309   }
 310
 311   fTPCsubEPres->Fill(evPlaneResTPC,cent);
 312
 313   Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),cent};
 314   fEPres->Fill(evPlaneRes);
 315
 316   // Pi0, eta and gamma weights
 317
 318   if(fIsMC && fMC && stack && cent>20 && cent<40){
 319    Int_t nParticles = stack->GetNtrack();
 320    for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
 321       TParticle* particle = stack->Particle(iParticle);
 322       int fPDG = particle->GetPdgCode();
 323       double pTMC = particle->Pt();
 324       double etaMC = particle->Eta();
 325       if(fabs(etaMC)>0.7)continue;
 326
 327       Bool_t MChijing = fMC->IsFromBGEvent(iParticle);
 328       int iHijing = 1;
 329       if(!MChijing)iHijing = 0;
 330
 331       if(fPDG==111)fPi0Weight->Fill(pTMC,iHijing);//pi0
 332       if(fPDG==221)fEtaWeight->Fill(pTMC,iHijing);//eta
 333
 334       Int_t idMother = particle->GetFirstMother();
 335       if (idMother>0){
 336         TParticle *mother = stack->Particle(idMother);
 337         int motherPDG = mother->GetPdgCode();
 338         if(fPDG==22 && motherPDG!=111 && motherPDG!=221)fGammaWeight->Fill(pTMC,iHijing);//gamma
 339       }
 340
 341    }
 342   }
 343
 344
 345   // Track loop
 346   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 347     AliESDtrack* track = fESD->GetTrack(iTracks);
 348     if (!track) {
 349       printf("ERROR: Could not receive track %d\n", iTracks);
 350       continue;
 351     }
 352
 353     if(TMath::Abs(track->Eta())>0.7) continue;
 354
 355     fTrackPtBefTrkCuts->Fill(track->Pt());
 356
 357     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 358
 359     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 360       if(track->GetKinkIndex(0) != 0) continue;
 361     }
 362
 363     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 364
 365     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 366
 367     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 368
 369     fTrackPtAftTrkCuts->Fill(track->Pt());
 370
 371     Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999., wclsE = -999., wEovP = -999.;//, m02= -999., m20= -999.;
 372
 373     pt = track->Pt();
 374     if(pt<2) continue;
 375     fTrkpt->Fill(pt);
 376
 377     Int_t clsId = track->GetEMCALcluster();
 378     if (clsId>0){
 379       AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
 380       if(cluster && cluster->IsEMCAL()){
 381         clsE = cluster->E();
 382 //      m20 = cluster->M20();
 383 //      m02 = cluster->M02();
 384       }
 385     }
 386
 387     p = track->P();
 388     phi = track->Phi();
 389     dEdx = track->GetTPCsignal();
 390     EovP = clsE/p;
 391     wEovP = wclsE/p;
 392     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 393     fdEdxBef->Fill(p,dEdx);
 394     fTPCnsigma->Fill(p,fTPCnSigma);
 395
 396     //Remove electron candidate from the event plane
 397     Float_t evPlaneCorrTPC = evPlaneTPC;
 398     if(dEdx>70 && dEdx<90){
 399       Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
 400       Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
 401       TVector2 newQVectorfortrack;
 402       newQVectorfortrack.Set(qX,qY);
 403       evPlaneCorrTPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
 404     }
 405
 406     Bool_t fFlagPhotonicElec = kFALSE;
 407     Bool_t fFlagPhotonicElecBCG = kFALSE;
 408
 409     SelectPhotonicElectron(iTracks,track, fFlagPhotonicElec, fFlagPhotonicElecBCG);
 410
 411     Double_t corr[10]={phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneCorrTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagPhotonicElecBCG};
 412     fCorr->Fill(corr);
 413
 414
 415     Int_t whichFirstMother = 0, whichSecondMother = 0, whichThirdMother = 0;
 416     Int_t whichPart = -99;
 417     Int_t partPDG = -99, motherPDG = -99, secondMotherPDG = -99, thirdMotherPDG = -99;
 418     Double_t partPt = -99. , motherPt = -99., secondMotherPt = -99.,thirdMotherPt = -99.;
 419     Bool_t MChijing;
 420
 421     if(fIsMC && fMC && stack){
 422       Int_t label = track->GetLabel();
 423       if(label>0){
 424         TParticle *particle = stack->Particle(label);
 425         if(particle){
 426             partPDG = particle->GetPdgCode();
 427             partPt = particle->Pt();
 428
 429             if (TMath::Abs(partPDG)==11) whichPart = 0; //electron
 430             if (partPDG==22) whichPart = 3; //gamma
 431             if (partPDG==111) whichPart = 2; //pi0
 432             if (partPDG==221) whichPart = 1; //eta
 433
 434             MChijing = fMC->IsFromBGEvent(label);
 435
 436             int iHijing = 1;
 437             if(!MChijing) iHijing = 0; // 0 if enhanced sample
 438
 439             Int_t idMother = particle->GetFirstMother();
 440             if (idMother>0){
 441               TParticle *mother = stack->Particle(idMother);
 442               motherPt = mother->Pt();
 443               motherPDG = mother->GetPdgCode();
 444
 445               if (motherPDG==22) whichFirstMother = 3; //gamma
 446               if (motherPDG==111) whichFirstMother = 2; //pi0
 447               if (motherPDG==221) whichFirstMother = 1; //eta
 448
 449               Int_t idSecondMother = particle->GetSecondMother();
 450               if (idSecondMother>0){
 451                 TParticle *secondMother = stack->Particle(idSecondMother);
 452                 secondMotherPt = secondMother->Pt();
 453                 secondMotherPDG = secondMother->GetPdgCode();
 454
 455                 if (secondMotherPDG==111) whichSecondMother = 2; //pi0
 456                 if (secondMotherPDG==221) whichSecondMother = 1; //eta
 457
 458                 Int_t idThirdMother = secondMother->GetFirstMother();
 459                 if (idThirdMother>0){
 460                   TParticle *thirdMother = stack->Particle(idThirdMother);
 461                   thirdMotherPt = thirdMother->Pt();
 462                   thirdMotherPDG = thirdMother->GetPdgCode();
 463
 464                   if (thirdMotherPDG==221) whichThirdMother = 1; //eta
 465                 }
 466               }
 467
 468               Double_t mc[15]={EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,whichPart,cent,pt,whichFirstMother,whichSecondMother,whichThirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt};
 469               fMCphotoElecPt->Fill(mc);
 470 //            if (motherPDG==22 || motherPDG==111 || motherPDG==221) fMCphotoElecPt->Fill(mc);// mother = gamma, pi0, eta
 471           }
 472         }
 473       }
 474     }
 475
 476
 477     if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
 478     Int_t pidpassed = 1;
 479
 480     //--- track accepted
 481     AliHFEpidObject hfetrack;
 482     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 483     hfetrack.SetRecTrack(track);
 484     hfetrack.SetPbPb();
 485     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 486
 487     Double_t corrV2[7]={phi,cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
 488     fChargPartV2->Fill(corrV2);
 489
 490     if(fTPCnSigma >= -0.5){
 491         Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
 492         feTPCV2->Fill(correctedV2);
 493     }
 494
 495     if(pidpassed==0) continue;
 496
 497     Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,evPlaneCorrTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
 498
 499     feV2->Fill(correctedV2);
 500
 501     fTrkEovPAft->Fill(pt,EovP);
 502     fdEdxAft->Fill(p,dEdx);
 503
 504   if(fFlagPhotonicElec){
 505     fphoteV2->Fill(correctedV2);
 506     fPhotoElecPt->Fill(pt);
 507   }
 508
 509   if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
 510
 511  }
 512  PostData(1, fOutputList);
 513 }
 514 //_________________________________________
 515 void AliAnalysisTaskElecV2::UserCreateOutputObjects()
 516 {
 517   //--- Check MC
 518   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
 519     fIsMC = kTRUE;
 520     printf("+++++ MC Data available");
 521   }
 522   //--------Initialize PID
 523   fPID->SetHasMCData(fIsMC);
 524
 525   if(!fPID->GetNumberOfPIDdetectors())
 526     {
 527       fPID->AddDetector("TPC", 0);
 528       fPID->AddDetector("EMCAL", 1);
 529     }
 530
 531   fPID->SortDetectors();
 532   fPIDqa = new AliHFEpidQAmanager();
 533   fPIDqa->Initialize(fPID);
 534
 535   //--------Initialize correction Framework and Cuts
 536   fCFM = new AliCFManager;
 537   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 538   fCFM->SetNStepParticle(kNcutSteps);
 539   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 540     fCFM->SetParticleCutsList(istep, NULL);
 541
 542   if(!fCuts){
 543     AliWarning("Cuts not available. Default cuts will be used");
 544     fCuts = new AliHFEcuts;
 545     fCuts->CreateStandardCuts();
 546   }
 547   fCuts->Initialize(fCFM);
 548
 549   //---------Output Tlist
 550   fOutputList = new TList();
 551   fOutputList->SetOwner();
 552   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 553
 554   fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
 555   fOutputList->Add(fNoEvents);
 556
 557   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
 558   fOutputList->Add(fTrkpt);
 559
 560   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
 561   fOutputList->Add(fTrackPtBefTrkCuts);
 562
 563   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
 564   fOutputList->Add(fTrackPtAftTrkCuts);
 565
 566   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
 567   fOutputList->Add(fTPCnsigma);
 568
 569   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
 570   fOutputList->Add(fTrkEovPBef);
 571
 572   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
 573   fOutputList->Add(fTrkEovPAft);
 574
 575   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
 576   fOutputList->Add(fdEdxBef);
 577
 578   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
 579   fOutputList->Add(fdEdxAft);
 580
 581   fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
 582   fOutputList->Add(fInvmassLS);
 583
 584   fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
 585   fOutputList->Add(fInvmassULS);
 586
 587   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
 588   fOutputList->Add(fOpeningAngleLS);
 589
 590   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
 591   fOutputList->Add(fOpeningAngleULS);
 592
 593   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
 594   fOutputList->Add(fPhotoElecPt);
 595
 596   fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
 597   fOutputList->Add(fSemiInclElecPt);
 598
 599   fCent = new TH1F("fCent","Centrality",100,0,100) ;
 600   fOutputList->Add(fCent);
 601
 602   fevPlaneV0A = new TH2F("fevPlaneV0A","V0A EP",100,0,TMath::Pi(),90,0,90);
 603   fOutputList->Add(fevPlaneV0A);
 604
 605   fevPlaneV0C = new TH2F("fevPlaneV0C","V0C EP",100,0,TMath::Pi(),90,0,90);
 606   fOutputList->Add(fevPlaneV0C);
 607
 608   fevPlaneTPC = new TH2F("fevPlaneTPC","TPC EP",100,0,TMath::Pi(),90,0,90);
 609   fOutputList->Add(fevPlaneTPC);
 610
 611   fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
 612   fOutputList->Add(fTPCsubEPres);
 613
 614   Int_t binsv1[4]={100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, cent
 615   Double_t xminv1[4]={-1,-1,-1,0};
 616   Double_t xmaxv1[4]={1,1,1,90};
 617   fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
 618   fOutputList->Add(fEPres);
 619
 620   //phi,fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C),fFlagPhotonicElec,fFlagPhotonicElecBCG,m20,m02
 621   Int_t binsv2[10]={100,100,90,100,100,100,100,100,3,3};
 622   Double_t xminv2[10]={0,-3.5,0,0,0,0,0,0,-1,-1};
 623   Double_t xmaxv2[10]={2*TMath::Pi(),3.5,90,50,3,TMath::Pi(),TMath::Pi(),TMath::Pi(),2,2};
 624   fCorr = new THnSparseD ("fCorr","Correlations",10,binsv2,xminv2,xmaxv2);
 625   fOutputList->Add(fCorr);
 626
 627   Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
 628   Double_t xminv3[5]={0,0,-1,-1,-1};
 629   Double_t xmaxv3[5]={90,50,1,1,1};
 630   feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
 631   fOutputList->Add(feV2);
 632
 633   Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 634   Double_t xminv4[5]={0,0,-1,-1,-1};
 635   Double_t xmaxv4[5]={90,50,1,1,1};
 636   fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
 637   fOutputList->Add(fphoteV2);
 638
 639   Int_t binsv5[7]={100,90,100,100,100,100,100}; // phi, cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 640   Double_t xminv5[7]={0,0,0,0,-1,-1,-1};
 641   Double_t xmaxv5[7]={2*TMath::Pi(),90,50,3,1,1,1};
 642   fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",7,binsv5,xminv5,xmaxv5);
 643   fOutputList->Add(fChargPartV2);
 644
 645   Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 646   Double_t xminv6[5]={0,0,-1,-1,-1};
 647   Double_t xmaxv6[5]={90,50,1,1,1};
 648   feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
 649   fOutputList->Add(feTPCV2);
 650
 651   //EovP,fTPCnSigma,partPt,fFlagPhotonicElec,fFlagPhotonicElecBCG,whichPart,cent,pt,firstMother,secondMother,thirdMother,iHijing,motherPt,secondMotherPt,thirdMotherPt
 652   Int_t binsv7[15]={100,100,100,3,3,5,90,100,5,5,5,3,100,100,100};
 653   Double_t xminv7[15]={0,-3.5,0,-1,-1,-1,0,0,-1,-1,-1,-1,0,0,0};
 654   Double_t xmaxv7[15]={3,3.5,50,2,2,4,90,50,4,4,4,2,50,50,50};
 655   fMCphotoElecPt = new THnSparseD ("fMCphotoElecPt", "pt distribution (MC)",15,binsv7,xminv7,xmaxv7);
 656   fOutputList->Add(fMCphotoElecPt);
 657
 658   fGammaWeight = new TH2F("fGammaWeight", "Gamma weight",100,0,50,3,-1,2);
 659   fOutputList->Add(fGammaWeight);
 660
 661   fPi0Weight = new TH2F("fPi0Weight", "Pi0 weight",100,0,50,3,-1,2);
 662   fOutputList->Add(fPi0Weight);
 663
 664   fEtaWeight = new TH2F("fEtaWeight", "Eta weight",100,0,50,3,-1,2);
 665   fOutputList->Add(fEtaWeight);
 666
 667   PostData(1,fOutputList);
 668 }
 669
 670 //________________________________________________________________________
 671 void AliAnalysisTaskElecV2::Terminate(Option_t *)
 672 {
 673   // Info("Terminate");
 674         AliAnalysisTaskSE::Terminate();
 675 }
 676
 677 //________________________________________________________________________
 678 Bool_t AliAnalysisTaskElecV2::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 679 {
 680   // Check single track cuts for a given cut step
 681   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 682   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 683   return kTRUE;
 684 }
 685 //_________________________________________
 686 void AliAnalysisTaskElecV2::SelectPhotonicElectron(Int_t iTracks,AliESDtrack *track,Bool_t &fFlagPhotonicElec, Bool_t &fFlagPhotonicElecBCG)
 687 {
 688   //Identify non-heavy flavour electrons using Invariant mass method
 689
 690   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 691   fTrackCuts->SetRequireTPCRefit(kTRUE);
 692   fTrackCuts->SetRequireITSRefit(kTRUE);
 693   fTrackCuts->SetEtaRange(-0.7,0.7);
 694   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 695   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
 696   fTrackCuts->SetMinNClustersTPC(100);
 697
 698   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 699
 700   Bool_t flagPhotonicElec = kFALSE;
 701   Bool_t flagPhotonicElecBCG = kFALSE;
 702
 703   for(Int_t jTracks = 0; jTracks<fESD->GetNumberOfTracks(); jTracks++){
 704
 705     if(jTracks==iTracks) continue;
 706
 707     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
 708     if (!trackAsso) {
 709       printf("ERROR: Could not receive track %d\n", jTracks);
 710       continue;
 711     }
 712
 713     Double_t dEdxAsso = -999., ptAsso=-999.;
 714     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 715     Double_t openingAngle = -999., mass=999., width = -999;
 716
 717     dEdxAsso = trackAsso->GetTPCsignal();
 718     ptAsso = trackAsso->Pt();
 719     Int_t chargeAsso = trackAsso->Charge();
 720     Int_t charge = track->Charge();
 721
 722     if(ptAsso <0.5) continue;
 723     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
 724     if(dEdxAsso <65 || dEdxAsso>100) continue;
 725
 726     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 727     if(charge>0) fPDGe1 = -11;
 728     if(chargeAsso>0) fPDGe2 = -11;
 729
 730     if(charge == chargeAsso) fFlagLS = kTRUE;
 731     if(charge != chargeAsso) fFlagULS = kTRUE;
 732
 733     AliKFParticle ge1(*track, fPDGe1);
 734     AliKFParticle ge2(*trackAsso, fPDGe2);
 735     AliKFParticle recg(ge1, ge2);
 736
 737     if(recg.GetNDF()<1) continue;
 738     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 739     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 740
 741     AliKFVertex primV(*pVtx);
 742     primV += recg;
 743     recg.SetProductionVertex(primV);
 744
 745     recg.SetMassConstraint(0,0.0001);
 746
 747     openingAngle = ge1.GetAngle(ge2);
 748     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
 749     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
 750
 751     if(openingAngle > fOpeningAngleCut) continue;
 752
 753     recg.GetMass(mass,width);
 754
 755     if(fFlagLS) fInvmassLS->Fill(mass);
 756     if(fFlagULS) fInvmassULS->Fill(mass);
 757
 758     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec) flagPhotonicElec = kTRUE;
 759     if(mass<fInvmassCut && fFlagLS && !flagPhotonicElecBCG) flagPhotonicElecBCG = kTRUE;
 760
 761   }
 762   fFlagPhotonicElec = flagPhotonicElec;
 763   fFlagPhotonicElecBCG = flagPhotonicElecBCG;
 764
 765 }
 766 //_________________________________________
 767 Double_t AliAnalysisTaskElecV2::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
 768 {
 769   //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
 770   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
 771   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
 772   Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
 773
 774   return cos2DeltaPhi;
 775 }
 776
 777 //_________________________________________
 778 Double_t AliAnalysisTaskElecV2::GetDeltaPhi(Double_t phiA,Double_t phiB) const
 779 {
 780   //Get phi-psi_EP
 781   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
 782   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
 783
 784   return dPhi;
 785 }