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 "AliAnalysisTaskHFECal.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 "AliCentrality.h"
  76
  77 ClassImp(AliAnalysisTaskHFECal)
  78 //________________________________________________________________________
  79 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name)
  80   : AliAnalysisTaskSE(name)
  81   ,fESD(0)
  82   ,fOutputList(0)
  83   ,fTrackCuts(0)
  84   ,fCuts(0)
  85   ,fIdentifiedAsOutInz(kFALSE)
  86   ,fPassTheEventCut(kFALSE)
  87   ,fRejectKinkMother(kFALSE)
  88   ,fVz(0.0)
  89   ,fCFM(0)
  90   ,fPID(0)
  91   ,fPIDqa(0)
  92   ,fOpeningAngleCut(0.1)
  93   ,fInvmassCut(0.01)
  94   ,fNoEvents(0)
  95   ,fEMCAccE(0)
  96   ,fTrkpt(0)
  97   ,fTrkEovPBef(0)
  98   ,fTrkEovPAft(0)
  99   ,fdEdxBef(0)
 100   ,fdEdxAft(0)
 101   ,fIncpT(0)
 102   ,fInvmassLS(0)
 103   ,fInvmassULS(0)
 104   ,fOpeningAngleLS(0)
 105   ,fOpeningAngleULS(0)
 106   ,fPhotoElecPt(0)
 107   ,fPhoElecPt(0)
 108   ,fTrackPtBefTrkCuts(0)
 109   ,fTrackPtAftTrkCuts(0)
 110   ,fTPCnsigma(0)
 111   ,fCent(0)
 112   ,fEleInfo(0)
 113 {
 114   //Named constructor
 115
 116   fPID = new AliHFEpid("hfePid");
 117   fTrackCuts = new AliESDtrackCuts();
 118
 119   // Define input and output slots here
 120   // Input slot #0 works with a TChain
 121   DefineInput(0, TChain::Class());
 122   // Output slot #0 id reserved by the base class for AOD
 123   // Output slot #1 writes into a TH1 container
 124   // DefineOutput(1, TH1I::Class());
 125   DefineOutput(1, TList::Class());
 126   //  DefineOutput(3, TTree::Class());
 127 }
 128
 129 //________________________________________________________________________
 130 AliAnalysisTaskHFECal::AliAnalysisTaskHFECal()
 131   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal")
 132   ,fESD(0)
 133   ,fOutputList(0)
 134   ,fTrackCuts(0)
 135   ,fCuts(0)
 136   ,fIdentifiedAsOutInz(kFALSE)
 137   ,fPassTheEventCut(kFALSE)
 138   ,fRejectKinkMother(kFALSE)
 139   ,fVz(0.0)
 140   ,fCFM(0)
 141   ,fPID(0)
 142   ,fPIDqa(0)
 143   ,fOpeningAngleCut(0.1)
 144   ,fInvmassCut(0.01)
 145   ,fNoEvents(0)
 146   ,fEMCAccE(0)
 147   ,fTrkpt(0)
 148   ,fTrkEovPBef(0)
 149   ,fTrkEovPAft(0)
 150   ,fdEdxBef(0)
 151   ,fdEdxAft(0)
 152   ,fIncpT(0)
 153   ,fInvmassLS(0)
 154   ,fInvmassULS(0)
 155   ,fOpeningAngleLS(0)
 156   ,fOpeningAngleULS(0)
 157   ,fPhotoElecPt(0)
 158   ,fPhoElecPt(0)
 159   ,fTrackPtBefTrkCuts(0)
 160   ,fTrackPtAftTrkCuts(0)
 161   ,fTPCnsigma(0)
 162   ,fCent(0)
 163   ,fEleInfo(0)
 164 {
 165         //Default constructor
 166         fPID = new AliHFEpid("hfePid");
 167
 168         fTrackCuts = new AliESDtrackCuts();
 169
 170         // Constructor
 171         // Define input and output slots here
 172         // Input slot #0 works with a TChain
 173         DefineInput(0, TChain::Class());
 174         // Output slot #0 id reserved by the base class for AOD
 175         // Output slot #1 writes into a TH1 container
 176         // DefineOutput(1, TH1I::Class());
 177         DefineOutput(1, TList::Class());
 178         //DefineOutput(3, TTree::Class());
 179 }
 180 //_________________________________________
 181
 182 AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal()
 183 {
 184   //Destructor
 185
 186   delete fOutputList;
 187   delete fPID;
 188   delete fCFM;
 189   delete fPIDqa;
 190   delete fTrackCuts;
 191 }
 192 //_________________________________________
 193
 194 void AliAnalysisTaskHFECal::UserExec(Option_t*)
 195 {
 196   //Main loop
 197   //Called for each event
 198
 199   // create pointer to event
 200   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 201   if (!fESD) {
 202     printf("ERROR: fESD not available\n");
 203     return;
 204   }
 205
 206   if(!fCuts){
 207     AliError("HFE cuts not available");
 208     return;
 209   }
 210
 211   if(!fPID->IsInitialized()){
 212     // Initialize PID with the given run number
 213     AliWarning("PID not initialised, get from Run no");
 214     fPID->InitializePID(fESD->GetRunNumber());
 215   }
 216
 217   fNoEvents->Fill(0);
 218
 219   Int_t fNOtrks =  fESD->GetNumberOfTracks();
 220   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 221
 222   Double_t pVtxZ = -999;
 223   pVtxZ = pVtx->GetZ();
 224
 225   if(TMath::Abs(pVtxZ)>10) return;
 226   fNoEvents->Fill(1);
 227
 228   if(fNOtrks<2) return;
 229   fNoEvents->Fill(2);
 230
 231   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 232   if(!pidResponse){
 233     AliDebug(1, "Using default PID Response");
 234     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 235   }
 236
 237   fPID->SetPIDResponse(pidResponse);
 238
 239   fCFM->SetRecEventInfo(fESD);
 240
 241   Float_t cent = -1.;
 242   AliCentrality *centrality = fESD->GetCentrality();
 243   cent = centrality->GetCentralityPercentile("V0M");
 244   fCent->Fill(cent);
 245
 246   if(cent>90.) return;
 247
 248  // Calorimeter info.
 249
 250   // make EMCAL array
 251   for(Int_t iCluster=0; iCluster<fESD->GetNumberOfCaloClusters(); iCluster++)
 252      {
 253       AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster);
 254       if(clust->IsEMCAL())
 255         {
 256          double clustE = clust->E();
 257          float  emcx[3]; // cluster pos
 258          clust->GetPosition(emcx);
 259          TVector3 clustpos(emcx[0],emcx[1],emcx[2]);
 260          double emcphi = clustpos.Phi();
 261          double emceta = clustpos.Eta();
 262          double calInfo[3];
 263          calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent;
 264          fEMCAccE->Fill(calInfo);
 265         }
 266    }
 267
 268   // Track loop
 269   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 270     AliESDtrack* track = fESD->GetTrack(iTracks);
 271     if (!track) {
 272       printf("ERROR: Could not receive track %d\n", iTracks);
 273       continue;
 274     }
 275
 276     if(TMath::Abs(track->Eta())>0.7) continue;
 277
 278     fTrackPtBefTrkCuts->Fill(track->Pt());
 279     // RecKine: ITSTPC cuts
 280     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 281
 282     //RecKink
 283     if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters
 284       if(track->GetKinkIndex(0) != 0) continue;
 285     }
 286
 287     // RecPrim
 288     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 289
 290     // HFEcuts: ITS layers cuts
 291     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 292
 293
 294     fTrackPtAftTrkCuts->Fill(track->Pt());
 295
 296     Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.;
 297
 298     // Track extrapolation
 299
 300     pt = track->Pt();
 301     if(pt<2.0)continue;
 302     Int_t charge = track->Charge();
 303     fTrkpt->Fill(pt);
 304     mom = track->P();
 305     phi = track->Phi();
 306     eta = track->Eta();
 307     dEdx = track->GetTPCsignal();
 308     fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 309
 310         double ncells = -1.0;
 311         double m20 = -1.0;
 312         double m02 = -1.0;
 313         double disp = -1.0;
 314         double rmatch = -1.0;
 315         double nmatch = -1.0;
 316
 317
 318     Int_t clsId = track->GetEMCALcluster();
 319     if (clsId>0){
 320       AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId);
 321       if(clust && clust->IsEMCAL()){
 322
 323                 double clustE = clust->E();
 324                 eop = clustE/fabs(mom);
 325                 //double clustT = clust->GetTOF();
 326                 ncells = clust->GetNCells();
 327                 m02 = clust->GetM02();
 328                 m20 = clust->GetM20();
 329                 disp = clust->GetDispersion();
 330                 double delphi = clust->GetTrackDx();
 331                 double deleta = clust->GetTrackDz();
 332                 rmatch = sqrt(pow(delphi,2)+pow(deleta,2));
 333                 nmatch = clust->GetNTracksMatched();
 334
 335                   double valdedx[14];
 336                   valdedx[0] = mom; valdedx[1] = pt; valdedx[2] = dEdx; valdedx[3] = phi; valdedx[4] = eta; valdedx[5] = fTPCnSigma;
 337                   valdedx[6] = eop; valdedx[7] = rmatch; valdedx[8] = ncells,  valdedx[9] = m02; valdedx[10] = m20; valdedx[11] = disp;
 338                   valdedx[12] = cent; valdedx[13] = charge;
 339                   fEleInfo->Fill(valdedx);
 340
 341
 342       }
 343     }
 344
 345     fdEdxBef->Fill(mom,dEdx);
 346     fTPCnsigma->Fill(mom,fTPCnSigma);
 347
 348
 349     // HFE cuts: TPC PID cleanup
 350     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 351
 352     if(fTPCnSigma >= 1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop);
 353     Int_t pidpassed = 1;
 354
 355
 356     //--- track accepted
 357     AliHFEpidObject hfetrack;
 358     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 359     hfetrack.SetRecTrack(track);
 360     hfetrack.SetPbPb();
 361     if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0;
 362
 363     if(pidpassed==0) continue;
 364
 365     fTrkEovPAft->Fill(pt,eop);
 366     fdEdxAft->Fill(mom,dEdx);
 367     fIncpT->Fill(cent,pt);
 368
 369     Bool_t fFlagPhotonicElec = kFALSE;
 370     SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec);
 371
 372     if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt);
 373  }
 374  PostData(1, fOutputList);
 375 }
 376 //_________________________________________
 377 void AliAnalysisTaskHFECal::UserCreateOutputObjects()
 378 {
 379   //--------Initialize PID
 380   fPID->SetHasMCData(kFALSE);
 381   if(!fPID->GetNumberOfPIDdetectors())
 382     {
 383       fPID->AddDetector("TPC", 0);
 384       fPID->AddDetector("EMCAL", 1);
 385     }
 386
 387   fPID->SortDetectors();
 388   fPIDqa = new AliHFEpidQAmanager();
 389   fPIDqa->Initialize(fPID);
 390
 391   //--------Initialize correction Framework and Cuts
 392   fCFM = new AliCFManager;
 393   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 394   fCFM->SetNStepParticle(kNcutSteps);
 395   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 396     fCFM->SetParticleCutsList(istep, NULL);
 397
 398   if(!fCuts){
 399     AliWarning("Cuts not available. Default cuts will be used");
 400     fCuts = new AliHFEcuts;
 401     fCuts->CreateStandardCuts();
 402   }
 403   fCuts->Initialize(fCFM);
 404
 405   //---------Output Tlist
 406   fOutputList = new TList();
 407   fOutputList->SetOwner();
 408   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 409
 410   fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ;
 411   fOutputList->Add(fNoEvents);
 412
 413   Int_t binsE[4] =    {250, 100,  1000, 100};
 414   Double_t xminE[4] = {1.0,  -1,   0.0,   0};
 415   Double_t xmaxE[4] = {3.5,   1, 100.0, 100};
 416   fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality",4,binsE,xminE,xmaxE);
 417   fOutputList->Add(fEMCAccE);
 418
 419   fTrkpt = new TH1F("fTrkpt","track pt",100,0,50);
 420   fOutputList->Add(fTrkpt);
 421
 422   fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50);
 423   fOutputList->Add(fTrackPtBefTrkCuts);
 424
 425   fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50);
 426   fOutputList->Add(fTrackPtAftTrkCuts);
 427
 428   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10);
 429   fOutputList->Add(fTPCnsigma);
 430
 431   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2);
 432   fOutputList->Add(fTrkEovPBef);
 433
 434   fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2);
 435   fOutputList->Add(fTrkEovPAft);
 436
 437   fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,150,0,150);
 438   fOutputList->Add(fdEdxBef);
 439
 440   fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,150,0,150);
 441   fOutputList->Add(fdEdxAft);
 442
 443   fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",100,0,100,100,0,50);
 444   fOutputList->Add(fIncpT);
 445
 446   fInvmassLS = new TH2F("fInvmassLS", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 100,0,100,500,0,0.5);
 447   fOutputList->Add(fInvmassLS);
 448
 449   fInvmassULS = new TH2F("fInvmassULS", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 100,0,100,500,0,0.5);
 450   fOutputList->Add(fInvmassULS);
 451
 452   fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
 453   fOutputList->Add(fOpeningAngleLS);
 454
 455   fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
 456   fOutputList->Add(fOpeningAngleULS);
 457
 458   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50);
 459   fOutputList->Add(fPhotoElecPt);
 460
 461   fPhoElecPt = new TH2F("fPhoElecPt", "Semi-inclusive electron pt",100,0,100,100,0,50);
 462   fOutputList->Add(fPhoElecPt);
 463
 464   fCent = new TH1F("fCent","Centrality",100,0,100) ;
 465   fOutputList->Add(fCent);
 466
 467   // Make common binning
 468   const Double_t kMinP = 0.;
 469   const Double_t kMaxP = 50.;
 470   const Double_t kTPCSigMim = 40.;
 471   const Double_t kTPCSigMax = 140.;
 472
 473   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta,  Sig,  e/p,  ,match, cell, M02, M20, Disp, Centrality, select)
 474   Int_t nBins[14] =  {  500,   500,        200,   60,    20,   600,  300, 100,   40,   200, 200, 200, 100,  3};
 475   Double_t min[14] = {kMinP, kMinP,  kTPCSigMim, 1.0,  -1.0,  -8.0,    0,   0,    0,   0.0, 0.0, 0.0,   0, -1.5};
 476   Double_t max[14] = {kMaxP, kMaxP,  kTPCSigMax, 4.0,   1.0,   4.0,  3.0, 0.1,   40,   2.0, 2.0, 2.0, 100,  1.5};
 477   fEleInfo = new THnSparseD("fEleInfo", "p [GeV/c]; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;M02;M20;Disp; Centrality; charge", 14, nBins, min, max);
 478   fOutputList->Add(fEleInfo);
 479
 480 //_________________________________________________________
 481
 482   PostData(1,fOutputList);
 483 }
 484
 485 //________________________________________________________________________
 486 void AliAnalysisTaskHFECal::Terminate(Option_t *)
 487 {
 488   // Info("Terminate");
 489         AliAnalysisTaskSE::Terminate();
 490 }
 491
 492 //________________________________________________________________________
 493 Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 494 {
 495   // Check single track cuts for a given cut step
 496   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 497   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 498   return kTRUE;
 499 }
 500 //_________________________________________
 501 void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec)
 502 {
 503   //Identify non-heavy flavour electrons using Invariant mass method
 504
 505   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 506   fTrackCuts->SetRequireTPCRefit(kTRUE);
 507   fTrackCuts->SetEtaRange(-0.7,0.7);
 508   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 509   fTrackCuts->SetMaxChi2PerClusterTPC(3.5);
 510   fTrackCuts->SetMinNClustersTPC(100);
 511
 512   const AliESDVertex *pVtx = fESD->GetPrimaryVertex();
 513
 514   Bool_t flagPhotonicElec = kFALSE;
 515
 516   for(Int_t jTracks = itrack+1; jTracks<fESD->GetNumberOfTracks(); jTracks++){
 517     AliESDtrack* trackAsso = fESD->GetTrack(jTracks);
 518     if (!trackAsso) {
 519       printf("ERROR: Could not receive track %d\n", jTracks);
 520       continue;
 521     }
 522
 523     Double_t dEdxAsso = -999., ptAsso=-999., openingAngle = -999.;
 524     Double_t mass=999., width = -999;
 525     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 526
 527     dEdxAsso = trackAsso->GetTPCsignal();
 528     ptAsso = trackAsso->Pt();
 529     Int_t chargeAsso = trackAsso->Charge();
 530     Int_t charge = track->Charge();
 531
 532     if(ptAsso <0.3) continue;
 533     if(!fTrackCuts->AcceptTrack(trackAsso)) continue;
 534     if(dEdxAsso <70 || dEdxAsso>100) continue; //11a pass1
 535
 536     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 537     if(charge>0) fPDGe1 = -11;
 538     if(chargeAsso>0) fPDGe2 = -11;
 539
 540     if(charge == chargeAsso) fFlagLS = kTRUE;
 541     if(charge != chargeAsso) fFlagULS = kTRUE;
 542
 543     AliKFParticle ge1(*track, fPDGe1);
 544     AliKFParticle ge2(*trackAsso, fPDGe2);
 545     AliKFParticle recg(ge1, ge2);
 546
 547     if(recg.GetNDF()<1) continue;
 548     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 549     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 550
 551     AliKFVertex primV(*pVtx);
 552     primV += recg;
 553     recg.SetProductionVertex(primV);
 554
 555     recg.SetMassConstraint(0,0.0001);
 556
 557     openingAngle = ge1.GetAngle(ge2);
 558     if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
 559     if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
 560
 561     if(openingAngle > fOpeningAngleCut) continue;
 562
 563     recg.GetMass(mass,width);
 564
 565     if(fFlagLS) fInvmassLS->Fill(cent,mass);
 566     if(fFlagULS) fInvmassULS->Fill(cent,mass);
 567
 568     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
 569       flagPhotonicElec = kTRUE;
 570     }
 571
 572   }
 573   fFlagPhotonicElec = flagPhotonicElec;
 574
 575 }
 576