PWGHF/hfe/AliAnalysisTaskFlowTPCEMCalQCSP.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
  17 ////////////////////////////////////////////////////////////////////////
  18 //                                                                    //
  19 //  Task for Heavy Flavour Electron Flow with TPC plus EMCal          //
  20 //  Non-Photonic Electron identified with Invariant mass              //
  21 //  analysis methos in function  SelectPhotonicElectron               //
  22 //                                                                    //
  23 //                                                                    //
  24 //  Author: Andrea Dubla (Utrecht University)                         //
  25 //                                                                                        //
  26 //                                                                    //
  27 ////////////////////////////////////////////////////////////////////////
  28
  29 #include "TChain.h"
  30 #include "TTree.h"
  31 #include "TH2F.h"
  32 #include "TMath.h"
  33 #include "TCanvas.h"
  34 #include "THnSparse.h"
  35 #include "TLorentzVector.h"
  36 #include "TString.h"
  37 #include "TFile.h"
  38 #include "AliAnalysisTask.h"
  39 #include "AliAnalysisManager.h"
  40 #include "AliESDEvent.h"
  41 #include "AliESDHandler.h"
  42 #include "AliAODEvent.h"
  43 #include "AliAODHandler.h"
  44 #include "AliAnalysisTaskFlowTPCEMCalQCSP.h"
  45 #include "TGeoGlobalMagField.h"
  46 #include "AliLog.h"
  47 #include "AliAnalysisTaskSE.h"
  48 #include "TRefArray.h"
  49 #include "TVector.h"
  50 #include "AliESDInputHandler.h"
  51 #include "AliESDpid.h"
  52 #include "AliAODInputHandler.h"
  53 #include "AliAODPid.h"
  54 #include "AliESDtrackCuts.h"
  55 #include "AliPhysicsSelection.h"
  56 #include "AliCentralitySelectionTask.h"
  57 #include "AliESDCaloCluster.h"
  58 #include "AliAODCaloCluster.h"
  59 #include "AliESDCaloTrigger.h"
  60 #include "AliEMCALRecoUtils.h"
  61 #include "AliEMCALGeometry.h"
  62 #include "AliGeomManager.h"
  63 #include "stdio.h"
  64 #include "TGeoManager.h"
  65 #include "iostream"
  66 #include "fstream"
  67 #include "AliEMCALTrack.h"
  68 //#include "AliEMCALTracker.h"
  69 #include "AliMagF.h"
  70 #include "AliKFParticle.h"
  71 #include "AliKFVertex.h"
  72 #include "AliPID.h"
  73 #include "AliPIDResponse.h"
  74 #include "AliHFEcontainer.h"
  75 #include "AliHFEcuts.h"
  76 #include "AliHFEpid.h"
  77 #include "AliHFEpidBase.h"
  78 #include "AliHFEpidQAmanager.h"
  79 #include "AliHFEtools.h"
  80 #include "AliCFContainer.h"
  81 #include "AliCFManager.h"
  82 #include "AliKFParticle.h"
  83 #include "AliKFVertex.h"
  84 #include "AliCentrality.h"
  85 #include "AliVEvent.h"
  86 #include "AliStack.h"
  87 #include "AliMCEvent.h"
  88 #include "TProfile.h"
  89 #include "AliFlowCandidateTrack.h"
  90 #include "AliFlowTrackCuts.h"
  91 #include "AliFlowEventSimple.h"
  92 #include "AliFlowCommonConstants.h"
  93 #include "AliFlowEvent.h"
  94 #include "TVector3.h"
  95 #include "TRandom2.h"
  96 #include "AliESDVZERO.h"
  97 #include "AliAODVZERO.h"
  98 #include "AliPID.h"
  99 #include "AliPIDResponse.h"
 100 #include "AliFlowTrack.h"
 101 #include "AliAnalysisTaskVnV0.h"
 102
 103
 104 class AliFlowTrackCuts;
 105
 106 using namespace std;
 107
 108 ClassImp(AliAnalysisTaskFlowTPCEMCalQCSP)
 109 //________________________________________________________________________
 110   AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP(const char *name)
 111   : AliAnalysisTaskSE(name)
 112 ,fDebug(0)
 113 ,fAOD(0)
 114 ,fGeom(0)
 115 ,fOutputList(0)
 116 ,fCuts(0)
 117 ,fIdentifiedAsOutInz(kFALSE)
 118 ,fPassTheEventCut(kFALSE)
 119 ,fCFM(0)
 120 ,fPID(0)
 121 ,fPIDqa(0)
 122 ,fCutsRP(0)     // track cuts for reference particles
 123 ,fNullCuts(0) // dummy cuts for flow event tracks
 124 ,fFlowEvent(0) //! flow events (one for each inv mass band)
 125 ,fkCentralityMethod(0)
 126 ,fCentrality(0)
 127 ,fCentralityMin(0)
 128 ,fCentralityMax(0)
 129 ,fInvmassCut(0)
 130 ,fTrigger(0)
 131 ,fPhi(0)
 132 ,fEta(0)
 133 ,fVZEROA(0)
 134 ,fVZEROC(0)
 135 ,fTPCM(0)
 136 ,fNoEvents(0)
 137 ,fTrkEovPBef(0)
 138 //,fdEdxBef(0)
 139 ,fInclusiveElecPt(0)
 140 ,fTPCnsigma(0)
 141 ,fTPCnsigmaAft(0)
 142 ,fCentralityPass(0)
 143 ,fCentralityNoPass(0)
 144 ,fInvmassLS1(0)
 145 ,fInvmassULS1(0)
 146 ,fPhotoElecPt(0)
 147 ,fSemiInclElecPt(0)
 148 ,fULSElecPt(0)
 149 ,fLSElecPt(0)
 150 ,fminTPC(-1)
 151 ,fmaxTPC(3)
 152 ,fminEovP(0.8)
 153 ,fmaxEovP(1.2)
 154 ,fminM20(0.03)
 155 ,fmaxM20(0.3)
 156 ,fminM02(0.03)
 157 ,fmaxM02(0.5)
 158 ,fDispersion(1)
 159 ,fMultCorAfterCuts(0)
 160 ,fMultvsCentr(0)
 161 ,fSubEventDPhiv2(0)
 162 ,EPVzA(0)
 163 ,EPVzC(0)
 164 ,EPTPC(0)
 165 ,fV2Phi(0)
 166 ,fSparseElectronHadron(0)
 167 ,fvertex(0)
 168 ,fMultCorBeforeCuts(0)
 169 ,fSideBandsFlow(kFALSE)
 170 ,fPhiminusPsi(kFALSE)
 171 ,fFlowEventCont(0) //! flow events (one for each inv mass band)
 172 ,fpurity(kFALSE)
 173 ,fSparseElectronpurity(0)
 174 {
 175   //Named constructor
 176
 177   fPID = new AliHFEpid("hfePid");
 178   // Define input and output slots here
 179   // Input slot #0 works with a TChain
 180   DefineInput(0, TChain::Class());
 181   // Output slot #0 id reserved by the base class for AOD
 182   // Output slot #1 writes into a TH1 container
 183   // DefineOutput(1, TH1I::Class());
 184   DefineOutput(1, TList::Class());
 185   DefineOutput(2, AliFlowEventSimple::Class());
 186 if(fSideBandsFlow){
 187     DefineOutput(3, AliFlowEventSimple::Class());
 188 }
 189  //  DefineOutput(3, TTree::Class());
 190 }
 191
 192 //________________________________________________________________________
 193 AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP()
 194   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskFlowTPCEMCalQCSP")
 195 ,fDebug(0)
 196 ,fAOD(0)
 197 ,fGeom(0)
 198 ,fOutputList(0)
 199 ,fCuts(0)
 200 ,fIdentifiedAsOutInz(kFALSE)
 201 ,fPassTheEventCut(kFALSE)
 202 ,fCFM(0)
 203 ,fPID(0)
 204 ,fPIDqa(0)
 205 ,fCutsRP(0)     // track cuts for reference particles
 206 ,fNullCuts(0) // dummy cuts for flow event tracks
 207 ,fFlowEvent(0) //! flow events (one for each inv mass band)
 208 ,fkCentralityMethod(0)
 209 ,fCentrality(0)
 210 ,fCentralityMin(0)
 211 ,fCentralityMax(0)
 212 ,fInvmassCut(0)
 213 ,fTrigger(0)
 214 ,fPhi(0)
 215 ,fEta(0)
 216 ,fVZEROA(0)
 217 ,fVZEROC(0)
 218 ,fTPCM(0)
 219 ,fNoEvents(0)
 220 ,fTrkEovPBef(0)
 221 //,fdEdxBef(0)
 222 ,fInclusiveElecPt(0)
 223 ,fTPCnsigma(0)
 224 ,fTPCnsigmaAft(0)
 225 ,fCentralityPass(0)
 226 ,fCentralityNoPass(0)
 227 ,fInvmassLS1(0)
 228 ,fInvmassULS1(0)
 229 ,fPhotoElecPt(0)
 230 ,fSemiInclElecPt(0)
 231 ,fULSElecPt(0)
 232 ,fLSElecPt(0)
 233 ,fminTPC(-1)
 234 ,fmaxTPC(3)
 235 ,fminEovP(0.8)
 236 ,fmaxEovP(1.2)
 237 ,fminM20(0.03)
 238 ,fmaxM20(0.3)
 239 ,fminM02(0.03)
 240 ,fmaxM02(0.5)
 241 ,fDispersion(1)
 242 ,fMultCorAfterCuts(0)
 243 ,fMultvsCentr(0)
 244 ,fSubEventDPhiv2(0)
 245 ,EPVzA(0)
 246 ,EPVzC(0)
 247 ,EPTPC(0)
 248 ,fV2Phi(0)
 249 ,fSparseElectronHadron(0)
 250 ,fvertex(0)
 251 ,fMultCorBeforeCuts(0)
 252 ,fSideBandsFlow(kFALSE)
 253 ,fPhiminusPsi(kFALSE)
 254 ,fFlowEventCont(0) //! flow events (one for each inv mass band)
 255 ,fpurity(kFALSE)
 256 ,fSparseElectronpurity(0)
 257 {
 258   //Default constructor
 259   fPID = new AliHFEpid("hfePid");
 260   // Constructor
 261   // Define input and output slots here
 262   // Input slot #0 works with a TChain
 263   DefineInput(0, TChain::Class());
 264   // Output slot #0 id reserved by the base class for AOD
 265   // Output slot #1 writes into a TH1 container
 266   // DefineOutput(1, TH1I::Class());
 267   DefineOutput(1, TList::Class());
 268   DefineOutput(2, AliFlowEventSimple::Class());
 269     //  DefineOutput(3, TTree::Class());
 270 if(fSideBandsFlow){
 271     DefineOutput(3, AliFlowEventSimple::Class());
 272 }
 273     //DefineOutput(3, TTree::Class());
 274 }
 275 //_________________________________________
 276
 277 AliAnalysisTaskFlowTPCEMCalQCSP::~AliAnalysisTaskFlowTPCEMCalQCSP()
 278 {
 279   //Destructor
 280
 281   delete fOutputList;
 282   delete fGeom;
 283   delete fPID;
 284   delete fCFM;
 285   delete fPIDqa;
 286   if (fOutputList) delete fOutputList;
 287   if (fFlowEvent) delete fFlowEvent;
 288   if (fFlowEventCont) delete fFlowEventCont;
 289
 290 }
 291 //_________________________________________
 292
 293 void AliAnalysisTaskFlowTPCEMCalQCSP::UserExec(Option_t*)
 294 {
 295   //Main loop
 296   //Called for each event
 297
 298   // create pointer to event
 299
 300   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 301   if (!fAOD)
 302   {
 303     printf("ERROR: fAOD not available\n");
 304     return;
 305   }
 306
 307   if(!fCuts)
 308   {
 309     AliError("HFE cuts not available");
 310     return;
 311   }
 312
 313   if(!fPID->IsInitialized())
 314   {
 315     // Initialize PID with the given run number
 316     AliWarning("PID not initialised, get from Run no");
 317     fPID->InitializePID(fAOD->GetRunNumber());
 318   }
 319
 320   // cout << "kTrigger   ==   " << fTrigger <<endl;
 321
 322     if(fTrigger==0){
 323 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral) ) return;
 324     }
 325     if(fTrigger==1){
 326 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
 327     }
 328     if(fTrigger==2){
 329 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA) ) return;
 330     }
 331     if(fTrigger==3){
 332 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) ) return;
 333     }
 334
 335
 336 //---------------CENTRALITY AND EVENT SELECTION-----------------------
 337    Int_t fNOtrks =  fAOD->GetNumberOfTracks();
 338     Float_t vtxz = -999;
 339     const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
 340     if (!trkVtx || trkVtx->GetNContributors()<=0)return;
 341     TString vtxTtl = trkVtx->GetTitle();
 342     if (!vtxTtl.Contains("VertexerTracks"))return;
 343     const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
 344     if (!spdVtx || spdVtx->GetNContributors()<=0)return;
 345     if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
 346     vtxz = trkVtx->GetZ();
 347     if(TMath::Abs(vtxz)>10)return;
 348     fvertex->Fill(vtxz);
 349
 350 // Event cut
 351     if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
 352     if(fNOtrks<2) return;
 353
 354     Bool_t pass = kFALSE; //to select centrality
 355     CheckCentrality(fAOD,pass);
 356     if(!pass)return;
 357
 358
 359   fNoEvents->Fill(0);
 360   PlotVZeroMultiplcities(fAOD);
 361
 362   SetNullCuts(fAOD);
 363   PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent);    //Calculate event plane Qvector and EP resolution for inclusive
 364
 365     if(fSideBandsFlow){
 366   PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEventCont);    //Calculate event plane Qvector and EP resolution for inclusive
 367     }
 368
 369   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 370   if(!pidResponse)
 371   {
 372     AliDebug(1, "Using default PID Response");
 373     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 374   }
 375
 376   fPID->SetPIDResponse(pidResponse);
 377   fCFM->SetRecEventInfo(fAOD);
 378
 379   // Look for kink mother
 380   Int_t numberofvertices = fAOD->GetNumberOfVertices();
 381   Double_t listofmotherkink[numberofvertices];
 382   Int_t numberofmotherkink = 0;
 383   for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
 384     AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
 385     if(!aodvertex) continue;
 386     if(aodvertex->GetType()==AliAODVertex::kKink) {
 387       AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
 388       if(!mother) continue;
 389       Int_t idmother = mother->GetID();
 390       listofmotherkink[numberofmotherkink] = idmother;
 391       //printf("ID %d\n",idmother);
 392       numberofmotherkink++;
 393     }
 394   }
 395
 396 //=============================================V0EP from Alex======================================================================
 397     Double_t qxEPa = 0, qyEPa = 0;
 398     Double_t qxEPc = 0, qyEPc = 0;
 399
 400     Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
 401     Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
 402
 403
 404     Double_t Qx2 = 0, Qy2 = 0;
 405
 406     for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
 407
 408         AliAODTrack* aodTrack = fAOD->GetTrack(iT);
 409
 410         if (!aodTrack)
 411             continue;
 412
 413         if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
 414             continue;
 415
 416         if (!aodTrack->TestFilterBit(128))
 417             continue;
 418
 419         Qx2 += TMath::Cos(2*aodTrack->Phi());
 420         Qy2 += TMath::Sin(2*aodTrack->Phi());
 421     }
 422
 423     Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
 424
 425     EPVzA->Fill(evPlAngV0A);
 426     EPVzC->Fill(evPlAngV0C);
 427     EPTPC->Fill(evPlAngTPC);
 428
 429     fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
 430     fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
 431     fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
 432 //====================================================================================================================
 433
 434
 435  AliAODTrack *track = NULL;
 436
 437 // Track loop
 438  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
 439  {
 440    track = fAOD->GetTrack(iTracks);
 441    if (!track)
 442    {
 443      printf("ERROR: Could not receive track %d\n", iTracks);
 444      continue;
 445    }
 446
 447    if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
 448 //----------hfe begin---------
 449    if(track->Eta()<-0.7 || track->Eta()>0.7)    continue;    //eta cuts on candidates
 450
 451    // RecKine: ITSTPC cuts
 452    if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 453
 454    // Reject kink mother
 455    Bool_t kinkmotherpass = kTRUE;
 456    for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
 457      if(track->GetID() == listofmotherkink[kinkmother]) {
 458        kinkmotherpass = kFALSE;
 459        continue;
 460      }
 461    }
 462    if(!kinkmotherpass) continue;
 463
 464    // RecPrim
 465  //  if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;  //deleted for DCA absence
 466    // HFEcuts: ITS layers cuts
 467    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 468    // HFE cuts: TPC PID cleanup
 469    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 470
 471    Double_t fClsE = -999, p = -999, fEovP=-999, pt = -999, fTPCnSigma=0;
 472    // Track extrapolation
 473    Int_t fClsId = track->GetEMCALcluster();
 474    if(fClsId < 0) continue;
 475    AliAODCaloCluster *cluster = fAOD->GetCaloCluster(fClsId);
 476    if(TMath::Abs(cluster->GetTrackDx()) > 0.05 || TMath::Abs(cluster->GetTrackDz()) > 0.05) continue;
 477
 478    pt = track->Pt();         //pt track after cuts
 479    if(pt<1) continue;
 480    fClsE = cluster->E();
 481    p = track->P();
 482   // dEdx = track->GetTPCsignal();
 483    fEovP = fClsE/p;
 484    fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 485    Double_t m20 =cluster->GetM20();
 486    Double_t m02 =cluster->GetM02();
 487    Double_t disp=cluster->GetDispersion();
 488    if(fTPCnSigma >= -1 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,fEovP);
 489    fTPCnsigma->Fill(p,fTPCnSigma);
 490 //   fdEdxBef->Fill(p,dEdx);
 491    Double_t eta = track->Eta();
 492    Double_t phi = track->Phi();
 493 //-----------------------Phiminupsi method to remove the contamination-----------------------------------------------
 494 //-----------------------fTPCnSigma < -3.5 hadrons will be selected from this region--------------------------
 495      Float_t dPhi_aeh = TVector2::Phi_0_2pi(phi - evPlAngV0A);
 496      if(dPhi_aeh > TMath::Pi()) dPhi_aeh = dPhi_aeh - TMath::Pi();
 497      Float_t dPhi_ceh = TVector2::Phi_0_2pi(phi - evPlAngV0C);
 498      if(dPhi_ceh > TMath::Pi()) dPhi_ceh = dPhi_ceh - TMath::Pi();
 499
 500      if(fPhiminusPsi){
 501          Double_t valueElh[8] = {
 502              pt,
 503              fEovP,
 504              fTPCnSigma,
 505              m20,
 506              m02,
 507              disp,
 508              dPhi_aeh,
 509              dPhi_ceh};
 510          fSparseElectronHadron->Fill(valueElh);
 511      }
 512 //----------------------------------------------------------------------------------------------------------
 513 //---------------------------From here usual electron selection---------------------------------------------
 514 //----------------------------------------------------------------------------------------------------------
 515 if(m20 < fminM20 || m20 > fmaxM20) continue;
 516 if(m02 < fminM02 || m02 > fmaxM02) continue;
 517 if(disp > fDispersion ) continue;
 518 //---------------------------------for purity---------------------------------------------------------------
 519      if(fpurity){
 520          Double_t valuepurity[3] = {
 521              pt,
 522              fEovP,
 523              fTPCnSigma};
 524          fSparseElectronpurity->Fill(valuepurity);
 525      }
 526 //----------------------------------------------------------------------------------------------------------
 527 //----------------------------------------------------------------------------------------------------------
 528 if(fTPCnSigma < fminTPC || fTPCnSigma > fmaxTPC) continue; //cuts on nsigma tpc and EoP
 529 //===============================Flow Event for Contamination=============================================
 530      if(fSideBandsFlow){
 531      if(fEovP>0 && fEovP<0.6){
 532          AliFlowTrack *sTrackCont = new AliFlowTrack();
 533          sTrackCont->Set(track);
 534          sTrackCont->SetID(track->GetID());
 535          sTrackCont->SetForRPSelection(kFALSE);
 536          sTrackCont->SetForPOISelection(kTRUE);
 537          sTrackCont->SetMass(2637);
 538          for(int iRPs=0; iRPs!=fFlowEventCont->NumberOfTracks(); ++iRPs)
 539          {
 540              //   cout << " no of rps " << iRPs << endl;
 541              AliFlowTrack *iRPCont = dynamic_cast<AliFlowTrack*>(fFlowEventCont->GetTrack( iRPs ));
 542              if (!iRPCont) continue;
 543              if (!iRPCont->InRPSelection()) continue;
 544              if( sTrackCont->GetID() == iRPCont->GetID())
 545              {
 546                  if(fDebug) printf(" was in RP set");
 547                  //       cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set" <<endl;
 548                  iRPCont->SetForRPSelection(kFALSE);
 549                  fFlowEventCont->SetNumberOfRPs(fFlowEventCont->GetNumberOfRPs() - 1);
 550              }
 551          } //end of for loop on RPs
 552          fFlowEventCont->InsertTrack(((AliFlowTrack*) sTrackCont));
 553      }
 554      }
 555 //==========================================================================================================
 556 //===============================From here eovP cut is used fro QC, SP and EPV0=============================
 557 if(fEovP < fminEovP || fEovP >fmaxEovP) continue;
 558 //==========================================================================================================
 559 //============================Event Plane Method with V0====================================================
 560      Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
 561      Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
 562      Double_t v2Phi[3] = {
 563          v2PhiV0A,
 564          v2PhiV0C,
 565          pt};
 566      fV2Phi->Fill(v2Phi);
 567 //=========================================================================================================
 568    fTPCnsigmaAft->Fill(p,fTPCnSigma);
 569    fInclusiveElecPt->Fill(pt);
 570    fPhi->Fill(phi);
 571    fEta->Fill(eta);
 572 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
 573    AliFlowTrack *sTrack = new AliFlowTrack();
 574      sTrack->Set(track);
 575      sTrack->SetID(track->GetID());
 576      sTrack->SetForRPSelection(kFALSE);
 577      sTrack->SetForPOISelection(kTRUE);
 578      sTrack->SetMass(263732);
 579      for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
 580      {
 581       //   cout << " no of rps " << iRPs << endl;
 582          AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
 583          if (!iRP) continue;
 584          if (!iRP->InRPSelection()) continue;
 585         if( sTrack->GetID() == iRP->GetID())
 586         {
 587           if(fDebug) printf(" was in RP set");
 588      //       cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set" <<endl;
 589             iRP->SetForRPSelection(kFALSE);
 590             fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
 591         }
 592       } //end of for loop on RPs
 593    fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
 594
 595 //----------------------Selection and Flow of Photonic Electrons-----------------------------
 596    Bool_t fFlagPhotonicElec = kFALSE;
 597    SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
 598    if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
 599    // Semi inclusive electron
 600    if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
 601
 602  }//end loop on track
 603
 604  PostData(1, fOutputList);
 605  PostData(2, fFlowEvent);
 606     if(fSideBandsFlow){
 607  PostData(3, fFlowEventCont);
 608     }
 609
 610  //----------hfe end---------
 611 }
 612 //_________________________________________
 613 void AliAnalysisTaskFlowTPCEMCalQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track, Bool_t &fFlagPhotonicElec)
 614 {
 615   //Identify non-heavy flavour electrons using Invariant mass method
 616
 617   Bool_t flagPhotonicElec = kFALSE;
 618
 619   for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
 620     AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
 621     if (!trackAsso) {
 622       printf("ERROR: Could not receive track %d\n", jTracks);
 623       continue;
 624     }
 625     //  if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
 626       if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
 627       if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)|| (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
 628
 629
 630     if(jTracks == itrack) continue;
 631     Double_t ptAsso=-999., nsigma=-999.0;
 632     Double_t mass=-999., width = -999;
 633     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 634
 635
 636     ptAsso = trackAsso->Pt();
 637     Short_t chargeAsso = trackAsso->Charge();
 638     Short_t charge = track->Charge();
 639     nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
 640
 641     if(trackAsso->GetTPCNcls() < 80) continue;
 642     if(nsigma < -3 || nsigma > 3) continue;
 643     if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
 644     if(ptAsso <0.3) continue;
 645
 646     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 647     if(charge>0) fPDGe1 = -11;
 648     if(chargeAsso>0) fPDGe2 = -11;
 649
 650     if(charge == chargeAsso) fFlagLS = kTRUE;
 651     if(charge != chargeAsso) fFlagULS = kTRUE;
 652
 653     AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
 654     AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
 655     AliKFParticle recg(ge1, ge2);
 656
 657     if(recg.GetNDF()<1) continue;
 658     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 659     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 660     recg.GetMass(mass,width);
 661
 662     if(fFlagLS) fInvmassLS1->Fill(mass);
 663     if(fFlagULS) fInvmassULS1->Fill(mass);
 664
 665        if(mass<fInvmassCut){
 666        if(fFlagULS){fULSElecPt->Fill(track->Pt());}
 667        if(fFlagLS){fLSElecPt->Fill(track->Pt());}
 668        }
 669
 670     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
 671       flagPhotonicElec = kTRUE;
 672     }
 673   }//track loop
 674   fFlagPhotonicElec = flagPhotonicElec;
 675 }
 676 //___________________________________________
 677 void AliAnalysisTaskFlowTPCEMCalQCSP::UserCreateOutputObjects()
 678 {
 679   //Create histograms
 680
 681   //----------hfe initialising begin---------
 682    fNullCuts = new AliFlowTrackCuts("null_cuts");
 683
 684    AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 685   cc->SetNbinsMult(10000);
 686   cc->SetMultMin(0);
 687   cc->SetMultMax(10000);
 688
 689   cc->SetNbinsPt(100);
 690   cc->SetPtMin(0);
 691   cc->SetPtMax(50);
 692
 693   cc->SetNbinsPhi(180);
 694   cc->SetPhiMin(0.0);
 695   cc->SetPhiMax(TMath::TwoPi());
 696
 697   cc->SetNbinsEta(200);
 698   cc->SetEtaMin(-7.0);
 699   cc->SetEtaMax(+7.0);
 700
 701   cc->SetNbinsQ(500);
 702   cc->SetQMin(0.0);
 703   cc->SetQMax(3.0);
 704
 705   //--------Initialize PID
 706   fPID->SetHasMCData(kFALSE);
 707   if(!fPID->GetNumberOfPIDdetectors())
 708   {
 709     fPID->AddDetector("TPC", 0);
 710     fPID->AddDetector("EMCAL", 1);
 711   }
 712
 713   fPID->SortDetectors();
 714   fPIDqa = new AliHFEpidQAmanager();
 715   fPIDqa->Initialize(fPID);
 716
 717   //--------Initialize correction Framework and Cuts
 718   fCFM = new AliCFManager;
 719   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 720   fCFM->SetNStepParticle(kNcutSteps);
 721   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 722     fCFM->SetParticleCutsList(istep, NULL);
 723
 724   if(!fCuts){
 725     AliWarning("Cuts not available. Default cuts will be used");
 726     fCuts = new AliHFEcuts;
 727     fCuts->CreateStandardCuts();
 728   }
 729
 730   fCuts->SetAOD();
 731   fCuts->Initialize(fCFM);
 732   //----------hfe initialising end--------
 733   //---------Output Tlist
 734   fOutputList = new TList();
 735   fOutputList->SetOwner();
 736   fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 737
 738   fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
 739   fOutputList->Add(fNoEvents);
 740
 741   fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",1000,0,50,200,-10,10);
 742   fOutputList->Add(fTPCnsigma);
 743
 744   fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",1000,0,50,200,-10,10);
 745   fOutputList->Add(fTPCnsigmaAft);
 746
 747   fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",1000,0,50,100,0,2);
 748   fOutputList->Add(fTrkEovPBef);
 749
 750 //    fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",1000,0,50,150,0,150);
 751  //     fOutputList->Add(fdEdxBef);
 752
 753   fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",1000,0,100);
 754   fOutputList->Add(fInclusiveElecPt);
 755
 756   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",1000,0,100);
 757   fOutputList->Add(fPhotoElecPt);
 758   //
 759   fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",1000,0,100);
 760   fOutputList->Add(fSemiInclElecPt);
 761
 762      fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",1000,0,100);
 763      fOutputList->Add(fULSElecPt);
 764
 765     fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",1000,0,100);
 766     fOutputList->Add(fLSElecPt);
 767
 768   fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
 769   fOutputList->Add(fInvmassLS1);
 770
 771   fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
 772   fOutputList->Add(fInvmassULS1);
 773
 774   fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
 775   fOutputList->Add(fCentralityPass);
 776
 777   fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
 778   fOutputList->Add(fCentralityNoPass);
 779
 780   fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
 781   fOutputList->Add(fPhi);
 782
 783   fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
 784   fOutputList->Add(fEta);
 785
 786   fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
 787   fOutputList->Add(fVZEROA);
 788
 789   fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
 790   fOutputList->Add(fVZEROC);
 791
 792   fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
 793   fOutputList->Add(fTPCM);
 794
 795   fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
 796   fOutputList->Add(fvertex);
 797
 798   fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
 799   fOutputList->Add(fMultCorBeforeCuts);
 800
 801   fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
 802   fOutputList->Add(fMultCorAfterCuts);
 803
 804   fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
 805   fOutputList->Add(fMultvsCentr);
 806
 807 //----------------------------------------------------------------------------
 808     EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
 809     fOutputList->Add(EPVzA);
 810     EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
 811     fOutputList->Add(EPVzC);
 812     EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
 813     fOutputList->Add(EPTPC);
 814 //----------------------------------------------------------------------------
 815     fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
 816     fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
 817     fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
 818     fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
 819     fOutputList->Add(fSubEventDPhiv2);
 820 //================================Event Plane with VZERO=====================
 821     const Int_t nPtBins = 10;
 822     Double_t binsPt[nPtBins+1] = {0, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
 823     // v2A, v2C, pt
 824     Int_t    bins[3] = {  50,  50, nPtBins};
 825     Double_t xmin[3] = { -1., -1.,   0};
 826     Double_t xmax[3] = {  1.,  1.,   8};
 827     fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
 828     // Set bin limits for axes which are not standard binned
 829     fV2Phi->SetBinEdges(2, binsPt);
 830     // set axes titles
 831     fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
 832     fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
 833     fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 834     fOutputList->Add(fV2Phi);
 835 //----------------------------------------------------------------------------
 836     if(fPhiminusPsi){
 837     Int_t binsvElectH[8]={ 600,  200, 200 ,100,  100,  100,   10,          10}; //pt, E/p,TPCnSigma,M20,M02,Disp Phi-psiV0A ,Phi-PsiV0C,eta (commented)
 838     Double_t xminvElectH[8]={0,    0, -10 ,  0,    0,    0,    0,           0};
 839     Double_t xmaxvElectH[8]={20,   2,  10 ,  2,    2,    2,  TMath::Pi(), TMath::Pi()};
 840     fSparseElectronHadron = new THnSparseD("ElectronHadron","ElectronHadron",8,binsvElectH,xminvElectH,xmaxvElectH);
 841     fSparseElectronHadron->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
 842     fSparseElectronHadron->GetAxis(1)->SetTitle("EovP");
 843     fSparseElectronHadron->GetAxis(2)->SetTitle("TPCnSigma");
 844     fSparseElectronHadron->GetAxis(3)->SetTitle("M20");
 845     fSparseElectronHadron->GetAxis(4)->SetTitle("M02");
 846     fSparseElectronHadron->GetAxis(5)->SetTitle("Disp");
 847     fSparseElectronHadron->GetAxis(6)->SetTitle("phiminuspsi V0A");
 848     fSparseElectronHadron->GetAxis(7)->SetTitle("phiminuspsi V0C");
 849     fOutputList->Add(fSparseElectronHadron);
 850     }
 851 //----------------------------------------------------------------------------
 852 //----------------------------------------------------------------------------
 853     if(fpurity){
 854     Int_t binsvpurity[3]={   600,200, 200}; //pt, E/p,TPCnSigma
 855     Double_t xminvpurity[3]={0,    0, -10};
 856     Double_t xmaxvpurity[3]={20,   2,  10};
 857     fSparseElectronpurity = new THnSparseD("Electronpurity","Electronpurity",3,binsvpurity,xminvpurity,xmaxvpurity);
 858     fSparseElectronpurity->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
 859     fSparseElectronpurity->GetAxis(1)->SetTitle("EovP");
 860     fSparseElectronpurity->GetAxis(2)->SetTitle("TPCnSigma");
 861     fOutputList->Add(fSparseElectronpurity);
 862     }
 863 //----------------------------------------------------------------------------
 864
 865   PostData(1,fOutputList);
 866  // create and post flowevent
 867   fFlowEvent = new AliFlowEvent(10000);
 868   PostData(2, fFlowEvent);
 869
 870     if(fSideBandsFlow){
 871     fFlowEventCont = new AliFlowEvent(10000);
 872     PostData(3, fFlowEventCont);
 873     }
 874  }
 875
 876 //________________________________________________________________________
 877 void AliAnalysisTaskFlowTPCEMCalQCSP::Terminate(Option_t *)
 878 {
 879   // Info("Terminate");
 880   AliAnalysisTaskSE::Terminate();
 881 }
 882 //_____________________________________________________________________________
 883 template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::PlotVZeroMultiplcities(const T* event) const
 884 {
 885   // QA multiplicity plots
 886   fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
 887   fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
 888 }
 889 //_____________________________________________________________________________
 890 template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::SetNullCuts(T* event)
 891 {
 892  //Set null cuts
 893     if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
 894     fCutsRP->SetEvent(event, MCEvent());
 895     fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
 896     fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
 897     fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
 898     fNullCuts->SetEvent(event, MCEvent());
 899 }
 900 //_____________________________________________________________________________
 901 void AliAnalysisTaskFlowTPCEMCalQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
 902 {
 903  //Prepare flow events
 904    FlowEv->ClearFast();
 905    FlowEv->Fill(fCutsRP, fNullCuts);
 906    FlowEv->SetReferenceMultiplicity(iMulti);
 907    FlowEv->DefineDeadZone(0, 0, 0, 0);
 908  //  FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
 909 }
 910 //_____________________________________________________________________________
 911 Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 912 {
 913   // Check single track cuts for a given cut step
 914   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 915   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 916   return kTRUE;
 917 }
 918 //_________________________________________
 919 void AliAnalysisTaskFlowTPCEMCalQCSP::CheckCentrality(AliAODEvent* event, Bool_t &centralitypass)
 920 {
 921   // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
 922   if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
 923   fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
 924 //  cout << "--------------Centrality evaluated-------------------------"<<endl;
 925   if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
 926   {
 927     fCentralityNoPass->Fill(fCentrality);
 928 //    cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
 929     centralitypass = kFALSE;
 930   }else
 931   {
 932 //    cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
 933     centralitypass = kTRUE;
 934   }
 935 //to remove the bias introduced by multeplicity outliers---------------------
 936     Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
 937     Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
 938
 939     if (TMath::Abs(centv0 - centTrk) > 5.0){
 940         centralitypass = kFALSE;
 941         fCentralityNoPass->Fill(fCentrality);
 942      }
 943     const Int_t nGoodTracks = event->GetNumberOfTracks();
 944
 945     Float_t multTPC(0.); // tpc mult estimate
 946     Float_t multGlob(0.); // global multiplicity
 947     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
 948         AliAODTrack* trackAOD = event->GetTrack(iTracks);
 949         if (!trackAOD) continue;
 950         if (!(trackAOD->TestFilterBit(1))) continue;
 951         if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70)  || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
 952         multTPC++;
 953     }
 954     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
 955         AliAODTrack* trackAOD = event->GetTrack(iTracks);
 956         if (!trackAOD) continue;
 957         if (!(trackAOD->TestFilterBit(16))) continue;
 958         if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
 959         Double_t b[2] = {-99., -99.};
 960         Double_t bCov[3] = {-99., -99., -99.};
 961         if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
 962         if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
 963         multGlob++;
 964     } //track loop
 965     //     printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
 966  //   if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))){  2010
 967     if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
 968         centralitypass = kFALSE;
 969         fCentralityNoPass->Fill(fCentrality);
 970     }//2011
 971     fMultCorBeforeCuts->Fill(multGlob, multTPC);
 972
 973     if(centralitypass){
 974     fCentralityPass->Fill(fCentrality);
 975     fMultCorAfterCuts->Fill(multGlob, multTPC);
 976     fMultvsCentr->Fill(fCentrality, multTPC);
 977     }
 978 }
 979 //_____________________________________________________________________________
 980 void AliAnalysisTaskFlowTPCEMCalQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
 981 {
 982   // Set a centrality range ]min, max] and define the method to use for centrality selection
 983   fCentralityMin = CentralityMin;
 984   fCentralityMax = CentralityMax;
 985   fkCentralityMethod = CentralityMethod;
 986 }
 987 //_____________________________________________________________________________
 988 void AliAnalysisTaskFlowTPCEMCalQCSP::SetIDCuts(Double_t minTPC, Double_t maxTPC, Double_t minEovP, Double_t maxEovP, Double_t minM20, Double_t maxM20, Double_t minM02, Double_t maxM02, Double_t Dispersion)
 989 {
 990     //Set ID cuts
 991     fminTPC = minTPC;
 992     fmaxTPC = maxTPC;
 993     fminEovP = minEovP;
 994     fmaxEovP = maxEovP;
 995     fminM20 = minM20;
 996     fmaxM20 = maxM20;
 997     fminM02 = minM02;
 998     fmaxM02 = maxM02;
 999     fDispersion = Dispersion;
1000 }
1001 //_____________________________________________________________________________
1002
1003