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 "AliPID.h"
  65 #include "AliPIDResponse.h"
  66 #include "AliHFEcontainer.h"
  67 #include "AliHFEcuts.h"
  68 #include "AliHFEpid.h"
  69 #include "AliHFEpidBase.h"
  70 #include "AliHFEpidQAmanager.h"
  71 #include "AliHFEtools.h"
  72 #include "AliCFContainer.h"
  73 #include "AliCFManager.h"
  74
  75 #include "AliEventplane.h"
  76 #include "AliCentrality.h"
  77
  78 ClassImp(AliAnalysisTaskElecV2)
  79 //________________________________________________________________________
  80 AliAnalysisTaskElecV2::AliAnalysisTaskElecV2(const char *name)
  81   : AliAnalysisTaskSE(name)
  82   ,fESD(0)
  83   ,fOutputList(0)
  84   ,fTrackCuts(0)
  85   ,fCuts(0)
  86   ,fIdentifiedAsOutInz(kFALSE)
  87   ,fPassTheEventCut(kFALSE)
  88   ,fRejectKinkMother(kFALSE)
  89   ,fVz(0.0)
  90   ,fCFM(0)
  91   ,fPID(0)
  92   ,fPIDqa(0)
  93   ,fOpeningAngleCut(0.1)
  94   ,fInvmassCut(0.01)
  95   ,fNoEvents(0)
  96   ,fTrkpt(0)
  97   ,fTrkEovPBef(0)
  98   ,fTrkEovPAft(0)
  99   ,fdEdxBef(0)
 100   ,fdEdxAft(0)
 101   ,fInvmassLS(0)
 102   ,fInvmassULS(0)
 103   ,fOpeningAngleLS(0)
 104   ,fOpeningAngleULS(0)
 105   ,fPhotoElecPt(0)
 106   ,fSemiInclElecPt(0)
 107   ,fTrackPtBefTrkCuts(0)
 108   ,fTrackPtAftTrkCuts(0)
 109   ,fTPCnsigma(0)
 110   ,fCent(0)
 111   ,fTPCsubEPres(0)
 112   ,fEPres(0)
 113   ,fCorr(0)
 114   ,feTPCV2(0)
 115   ,feV2(0)
 116   ,fphoteV2(0)
 117   ,fChargPartV2(0)
 118 {
 119   //Named constructor
 120
 121   fPID = new AliHFEpid("hfePid");
 122   fTrackCuts = new AliESDtrackCuts();
 123
 124   // Define input and output slots here
 125   // Input slot #0 works with a TChain
 126   DefineInput(0, TChain::Class());
 127   // Output slot #0 id reserved by the base class for AOD
 128   // Output slot #1 writes into a TH1 container
 129   // DefineOutput(1, TH1I::Class());
 130   DefineOutput(1, TList::Class());
 131   //  DefineOutput(3, TTree::Class());
 132 }
 133
 134 //________________________________________________________________________
 135 AliAnalysisTaskElecV2::AliAnalysisTaskElecV2()
 136   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElecHadCorrel")
 137   ,fESD(0)
 138   ,fOutputList(0)
 139   ,fTrackCuts(0)
 140   ,fCuts(0)
 141   ,fIdentifiedAsOutInz(kFALSE)
 142   ,fPassTheEventCut(kFALSE)
 143   ,fRejectKinkMother(kFALSE)
 144   ,fVz(0.0)
 145   ,fCFM(0)
 146   ,fPID(0)
 147   ,fPIDqa(0)
 148   ,fOpeningAngleCut(0.1)
 149   ,fInvmassCut(0.01)
 150   ,fNoEvents(0)
 151   ,fTrkpt(0)
 152   ,fTrkEovPBef(0)
 153   ,fTrkEovPAft(0)
 154   ,fdEdxBef(0)
 155   ,fdEdxAft(0)
 156   ,fInvmassLS(0)
 157   ,fInvmassULS(0)
 158   ,fOpeningAngleLS(0)
 159   ,fOpeningAngleULS(0)
 160   ,fPhotoElecPt(0)
 161   ,fSemiInclElecPt(0)
 162   ,fTrackPtBefTrkCuts(0)
 163   ,fTrackPtAftTrkCuts(0)
 164   ,fTPCnsigma(0)
 165   ,fCent(0)
 166   ,fTPCsubEPres(0)
 167   ,fEPres(0)
 168   ,fCorr(0)
 169   ,feTPCV2(0)
 170   ,feV2(0)
 171   ,fphoteV2(0)
 172   ,fChargPartV2(0)
 173 {
 174         //Default constructor
 175         fPID = new AliHFEpid("hfePid");
 176
 177         fTrackCuts = new AliESDtrackCuts();
 178
 179         // Constructor
 180         // Define input and output slots here
 181         // Input slot #0 works with a TChain
 182         DefineInput(0, TChain::Class());
 183         // Output slot #0 id reserved by the base class for AOD
 184         // Output slot #1 writes into a TH1 container
 185         // DefineOutput(1, TH1I::Class());
 186         DefineOutput(1, TList::Class());
 187         //DefineOutput(3, TTree::Class());
 188 }
 189 //_________________________________________
 190
 191 AliAnalysisTaskElecV2::~AliAnalysisTaskElecV2()
 192 {
 193   //Destructor
 194
 195   delete fOutputList;
 196   delete fPID;
 197   delete fCFM;
 198   delete fPIDqa;
 199   delete fTrackCuts;
 200 }
 201 //_________________________________________
 202
 203 void AliAnalysisTaskElecV2::UserExec(Option_t*)
 204 {
 205   //Main loop
 206   //Called for each event
 207
 208   // create pointer to event
 209   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 210   if (!fESD) {
 211     printf("ERROR: fESD not available\n");
 212     return;
 213   }
 214
 215   if(!fCuts){
 216     AliError("HFE cuts not available");
 217     return;
 218   }
 219
 220   if(!fPID->IsInitialized()){
 221     // Initialize PID with the given run number
 222     AliWarning("PID not initialised, get from Run no");
 223     fPID->InitializePID(fESD->GetRunNumber());
 224   }
 225
 226
 227   Int_t fNOtrks =  fESD->GetNumberOfTracks();
 228   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 229
 230   Double_t pVtxZ = -999;
 231   pVtxZ = pVtx->GetZ();
 232
 233   if(TMath::Abs(pVtxZ)>10) return;
 234   fNoEvents->Fill(0);
 235
 236   if(fNOtrks<2) return;
 237
 238   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 239   if(!pidResponse){
 240     AliDebug(1, "Using default PID Response");
 241     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 242   }
 243
 244   fPID->SetPIDResponse(pidResponse);
 245
 246   fCFM->SetRecEventInfo(fESD);
 247
 248   Float_t cent = -1.;
 249   AliCentrality *centrality = fESD->GetCentrality();
 250   cent = centrality->GetCentralityPercentile("V0M");
 251   fCent->Fill(cent);
 252
 253   if(cent>90.) return;
 254
 255   //Event planes
 256
 257   Double_t evPlaneV0A = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0A",fESD,2));
 258   if(evPlaneV0A > TMath::Pi()) evPlaneV0A = evPlaneV0A - TMath::Pi();
 259
 260   Double_t evPlaneV0C = TVector2::Phi_0_2pi(fESD->GetEventplane()->GetEventplane("V0C",fESD,2));
 261   if(evPlaneV0C > TMath::Pi()) evPlaneV0C = evPlaneV0C - TMath::Pi();
 262
 263   AliEventplane* esdTPCep = fESD->GetEventplane();
 264   TVector2 *standardQ = esdTPCep->GetQVector();
 265   Double_t qx = -999., qy = -999.;
 266   if(standardQ)
 267   {
 268           qx = standardQ->X();
 269           qy = standardQ->Y();
 270   }
 271   TVector2 qVectorfortrack;
 272   qVectorfortrack.Set(qx,qy);
 273   Float_t evPlaneTPC = TVector2::Phi_0_2pi(qVectorfortrack.Phi())/2.;
 274
 275   TVector2 *qsub1a = esdTPCep->GetQsub1();
 276   TVector2 *qsub2a = esdTPCep->GetQsub2();
 277   Double_t evPlaneResTPC = -999.;
 278   if(qsub1a && qsub2a)
 279   {
 280           evPlaneResTPC = TMath::Cos(2.*TVector2::Phi_0_2pi(qsub1a->Phi()/2.- qsub2a->Phi()/2.));
 281   }
 282
 283   fTPCsubEPres->Fill(evPlaneResTPC,cent);
 284
 285   Double_t evPlaneRes[4]={GetCos2DeltaPhi(evPlaneV0A,evPlaneV0C),GetCos2DeltaPhi(evPlaneV0A,evPlaneTPC),GetCos2DeltaPhi(evPlaneV0C,evPlaneTPC),cent};
 286   fEPres->Fill(evPlaneRes);
 287
 288   // Track loop
 289   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 290     AliESDtrack* track = fESD->GetTrack(iTracks);
 291     if (!track) {
 292       printf("ERROR: Could not receive track %d\n", iTracks);
 293       continue;
 294     }
 295
 296     if(TMath::Abs(track->Eta())>0.7) continue;
 297
 298     fTrackPtBefTrkCuts->Fill(track->Pt());
 299     // RecKine: ITSTPC cuts
 300     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 301
 302     //RecKink
 303     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 304       if(track->GetKinkIndex(0) != 0) continue;
 305     }
 306
 307     // RecPrim
 308     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 309
 310     // HFEcuts: ITS layers cuts
 311     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 312
 313     // HFE cuts: TPC PID cleanup
 314     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 315
 316     fTrackPtAftTrkCuts->Fill(track->Pt());
 317
 318     Double_t clsE = -999., p = -999., EovP=-999., pt = -999., dEdx=-999., fTPCnSigma=0, phi=-999.;
 319
 320     // Track extrapolation
 321
 322     pt = track->Pt();
 323     fTrkpt->Fill(pt);
 324
 325     Int_t clsId = track->GetEMCALcluster();
 326     if (clsId>0){
 327       AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsId);
 328       if(cluster && cluster->IsEMCAL()){
 329         clsE = cluster->E();
 330       }
 331     }
 332
 333     p = track->P();
 334     phi = track->Phi();
 335     dEdx = track->GetTPCsignal();
 336     EovP = clsE/p;
 337     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 338     fdEdxBef->Fill(p,dEdx);
 339     fTPCnsigma->Fill(p,fTPCnSigma);
 340
 341     Double_t corr[7]={fTPCnSigma,cent,pt,EovP,GetDeltaPhi(phi,evPlaneTPC),GetDeltaPhi(phi,evPlaneV0A),GetDeltaPhi(phi,evPlaneV0C)};
 342     fCorr->Fill(corr);
 343
 344     if(fTPCnSigma >= 1.5 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,EovP);
 345     Int_t pidpassed = 1;
 346
 347     //--- track accepted
 348     AliHFEpidObject hfetrack;
 349     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 350     hfetrack.SetRecTrack(track);
 351     hfetrack.SetPbPb();
 352     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 353
 354     Double_t corrV2[6]={cent,pt,EovP,GetCos2DeltaPhi(phi,evPlaneTPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
 355     fChargPartV2->Fill(corrV2);
 356
 357     if(fTPCnSigma >= -0.5){
 358         Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
 359         Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
 360         TVector2 newQVectorfortrack;
 361         newQVectorfortrack.Set(qX,qY);
 362         Double_t corrV2TPC = -999.;
 363         corrV2TPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
 364
 365         Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,corrV2TPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
 366
 367         feTPCV2->Fill(correctedV2);
 368     }
 369
 370     if(pidpassed==0) continue;
 371
 372     Double_t qX = standardQ->X() - esdTPCep->GetQContributionX(track);
 373     Double_t qY = standardQ->Y() - esdTPCep->GetQContributionY(track);
 374     TVector2 newQVectorfortrack;
 375     newQVectorfortrack.Set(qX,qY);
 376     Double_t corrV2TPC = -999.;
 377     corrV2TPC = TVector2::Phi_0_2pi(newQVectorfortrack.Phi())/2;
 378
 379     Double_t correctedV2[5]={cent,pt,GetCos2DeltaPhi(phi,corrV2TPC),GetCos2DeltaPhi(phi,evPlaneV0A),GetCos2DeltaPhi(phi,evPlaneV0C)};
 380
 381     feV2->Fill(correctedV2);
 382
 383     fTrkEovPAft->Fill(pt,EovP);
 384     fdEdxAft->Fill(p,dEdx);
 385
 386     Bool_t fFlagPhotonicElec = kFALSE;
 387     SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
 388
 389     if(fFlagPhotonicElec){
 390       fphoteV2->Fill(correctedV2);
 391       fPhotoElecPt->Fill(pt);
 392     }
 393
 394     if(!fFlagPhotonicElec) fSemiInclElecPt->Fill(pt);
 395  }
 396  PostData(1, fOutputList);
 397 }
 398 //_________________________________________
 399 void AliAnalysisTaskElecV2::UserCreateOutputObjects()
 400 {
 401   //--------Initialize PID
 402   fPID->SetHasMCData(kFALSE);
 403   if(!fPID->GetNumberOfPIDdetectors())
 404     {
 405       fPID->AddDetector("TPC", 0);
 406       fPID->AddDetector("EMCAL", 1);
 407     }
 408
 409   fPID->SortDetectors();
 410   fPIDqa = new AliHFEpidQAmanager();
 411   fPIDqa->Initialize(fPID);
 412
 413   //--------Initialize correction Framework and Cuts
 414   fCFM = new AliCFManager;
 415   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 416   fCFM->SetNStepParticle(kNcutSteps);
 417   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 418     fCFM->SetParticleCutsList(istep, NULL);
 419
 420   if(!fCuts){
 421     AliWarning("Cuts not available. Default cuts will be used");
 422     fCuts = new AliHFEcuts;
 423     fCuts->CreateStandardCuts();
 424   }
 425   fCuts->Initialize(fCFM);
 426
 427   //---------Output Tlist
 428   fOutputList = new TList();
 429   fOutputList->SetOwner();
 430   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 431
 432   fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
 433   fOutputList->Add(fNoEvents);
 434
 435   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
 436   fOutputList->Add(fTrkpt);
 437
 438   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
 439   fOutputList->Add(fTrackPtBefTrkCuts);
 440
 441   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
 442   fOutputList->Add(fTrackPtAftTrkCuts);
 443
 444   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
 445   fOutputList->Add(fTPCnsigma);
 446
 447   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
 448   fOutputList->Add(fTrkEovPBef);
 449
 450   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
 451   fOutputList->Add(fTrkEovPAft);
 452
 453   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
 454   fOutputList->Add(fdEdxBef);
 455
 456   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
 457   fOutputList->Add(fdEdxAft);
 458
 459   fInvmassLS = new TH1F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
 460   fOutputList->Add(fInvmassLS);
 461
 462   fInvmassULS = new TH1F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 500,0,0.5);
 463   fOutputList->Add(fInvmassULS);
 464
 465   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
 466   fOutputList->Add(fOpeningAngleLS);
 467
 468   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
 469   fOutputList->Add(fOpeningAngleULS);
 470
 471   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
 472   fOutputList->Add(fPhotoElecPt);
 473
 474   fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",100,0,50);
 475   fOutputList->Add(fSemiInclElecPt);
 476
 477   fCent = new TH1F("fCent","Centrality",100,0,100) ;
 478   fOutputList->Add(fCent);
 479
 480   fTPCsubEPres = new TH2F("fTPCsubEPres","TPC subevent plane resolution",100,-1,1,90,0,90);
 481   fOutputList->Add(fTPCsubEPres);
 482
 483   Int_t binsv1[4]={100,100,100,90}; // V0A-V0C, V0A-TPC, V0C-TPC, cent
 484   Double_t xminv1[4]={-1,-1,-1,0};
 485   Double_t xmaxv1[4]={1,1,1,90};
 486   fEPres = new THnSparseD ("fEPres","EP resolution",4,binsv1,xminv1,xmaxv1);
 487   fOutputList->Add(fEPres);
 488
 489   Int_t binsv2[7]={100,90,100,100,100,100,100}; // fTPCnSigma,cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 490   Double_t xminv2[7]={-3.5,0,0,0,0,0,0};
 491   Double_t xmaxv2[7]={3.5,90,50,3,TMath::Pi(),TMath::Pi(),TMath::Pi()};
 492   fCorr = new THnSparseD ("fCorr","Correlations",7,binsv2,xminv2,xmaxv2);
 493   fOutputList->Add(fCorr);
 494
 495   Int_t binsv3[5]={90,100,100,100,100}; // cent, pt, TPCcos2DeltaPhi, V0Acos2DeltaPhi, V0Ccos2DeltaPhi
 496   Double_t xminv3[5]={0,0,-1,-1,-1};
 497   Double_t xmaxv3[5]={90,50,1,1,1};
 498   feV2 = new THnSparseD ("feV2","inclusive electron v2",5,binsv3,xminv3,xmaxv3);
 499   fOutputList->Add(feV2);
 500
 501   Int_t binsv4[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 502   Double_t xminv4[5]={0,0,-1,-1,-1};
 503   Double_t xmaxv4[5]={90,50,1,1,1};
 504   fphoteV2 = new THnSparseD ("fphoteV2","photonic electron v2",5,binsv4,xminv4,xmaxv4);
 505   fOutputList->Add(fphoteV2);
 506
 507   Int_t binsv5[6]={90,100,100,100,100,100}; // cent, pt, EovP, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 508   Double_t xminv5[6]={0,0,0,-1,-1,-1};
 509   Double_t xmaxv5[6]={90,50,3,1,1,1};
 510   fChargPartV2 = new THnSparseD ("fChargPartV2","Charged particle v2",6,binsv5,xminv5,xmaxv5);
 511   fOutputList->Add(fChargPartV2);
 512
 513   Int_t binsv6[5]={90,100,100,100,100}; // cent, pt, TPCdeltaPhi, V0AdeltaPhi, V0CdeltaPhi
 514   Double_t xminv6[5]={0,0,-1,-1,-1};
 515   Double_t xmaxv6[5]={90,50,1,1,1};
 516   feTPCV2 = new THnSparseD ("feTPCV2","inclusive electron v2 (TPC)",5,binsv6,xminv6,xmaxv6);
 517   fOutputList->Add(feTPCV2);
 518
 519   PostData(1,fOutputList);
 520 }
 521
 522 //________________________________________________________________________
 523 void AliAnalysisTaskElecV2::Terminate(Option_t *)
 524 {
 525   // Info("Terminate");
 526         AliAnalysisTaskSE::Terminate();
 527 }
 528
 529 //________________________________________________________________________
 530 Bool_t AliAnalysisTaskElecV2::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 531 {
 532   // Check single track cuts for a given cut step
 533   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 534   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 535   return kTRUE;
 536 }
 537 //_________________________________________
 538 void AliAnalysisTaskElecV2::SelectPhotonicElectron(Int_t itrack, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
 539 {
 540   //Identify non-heavy flavour electrons using Invariant mass method
 541
 542   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 543   fTrackCuts->SetRequireTPCRefit(kTRUE);
 544   fTrackCuts->SetEtaRange(-0.7,0.7);
 545   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 546   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
 547   fTrackCuts->SetMinNClustersTPC(100);
 548
 549   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 550
 551   Bool_t flagPhotonicElec = kFALSE;
 552
 553   for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
 554     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
 555     if (!trackAsso) {
 556       printf("ERROR: Could not receive track %d\n", jTracks);
 557       continue;
 558     }
 559
 560     Double_t dEdxAsso = -999., ptAsso=-999., openingAngle = -999.;
 561     Double_t mass=999., width = -999;
 562     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 563
 564     dEdxAsso = trackAsso->GetTPCsignal();
 565     ptAsso = trackAsso->Pt();
 566     Int_t chargeAsso = trackAsso->Charge();
 567     Int_t charge = track->Charge();
 568
 569     if(ptAsso <0.3) continue;
 570     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
 571     if(dEdxAsso <70 || dEdxAsso>100) continue; //11a pass1
 572
 573     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 574     if(charge>0) fPDGe1 = -11;
 575     if(chargeAsso>0) fPDGe2 = -11;
 576
 577     if(charge == chargeAsso) fFlagLS = kTRUE;
 578     if(charge != chargeAsso) fFlagULS = kTRUE;
 579
 580     AliKFParticle ge1(*track, fPDGe1);
 581     AliKFParticle ge2(*trackAsso, fPDGe2);
 582     AliKFParticle recg(ge1, ge2);
 583
 584     if(recg.GetNDF()<1) continue;
 585     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 586     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 587
 588     AliKFVertex primV(*pVtx);
 589     primV += recg;
 590     recg.SetProductionVertex(primV);
 591
 592     recg.SetMassConstraint(0,0.0001);
 593
 594     openingAngle = ge1.GetAngle(ge2);
 595     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
 596     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
 597
 598     if(openingAngle > fOpeningAngleCut) continue;
 599
 600     recg.GetMass(mass,width);
 601
 602     if(fFlagLS) fInvmassLS->Fill(mass);
 603     if(fFlagULS) fInvmassULS->Fill(mass);
 604
 605     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
 606       flagPhotonicElec = kTRUE;
 607     }
 608
 609   }
 610   fFlagPhotonicElec = flagPhotonicElec;
 611
 612 }
 613 //_________________________________________
 614 Double_t AliAnalysisTaskElecV2::GetCos2DeltaPhi(Double_t phiA,Double_t phiB) const
 615 {
 616   //Get cos[2(phi-psi_EP)] or cos[2(psi_subEP1 - psi_subEP2)]
 617   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
 618   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
 619   Double_t cos2DeltaPhi = TMath::Cos(2*dPhi);
 620
 621   return cos2DeltaPhi;
 622 }
 623
 624 //_________________________________________
 625 Double_t AliAnalysisTaskElecV2::GetDeltaPhi(Double_t phiA,Double_t phiB) const
 626 {
 627   //Get phi-psi_EP
 628   Double_t dPhi = TVector2::Phi_0_2pi(phiA - phiB);
 629   if(dPhi > TMath::Pi()) dPhi = dPhi - TMath::Pi();
 630
 631   return dPhi;
 632 }