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 #include "AliSelectNonHFE.h"
 103
 104
 105 class AliFlowTrackCuts;
 106
 107 using namespace std;
 108
 109 ClassImp(AliAnalysisTaskFlowTPCEMCalQCSP)
 110 //________________________________________________________________________
 111 AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP(const char *name)
 112 : AliAnalysisTaskSE(name)
 113 ,fDebug(0)
 114 ,fAOD(0)
 115 ,fGeom(0)
 116 ,fOutputList(0)
 117 ,fCuts(0)
 118 ,fIdentifiedAsOutInz(kFALSE)
 119 ,fPassTheEventCut(kFALSE)
 120 ,fCFM(0)
 121 ,fPID(0)
 122 ,fPIDqa(0)
 123 ,fCutsRP(0)     // track cuts for reference particles
 124 ,fNullCuts(0) // dummy cuts for flow event tracks
 125 ,fFlowEvent(0) //! flow events (one for each inv mass band)
 126 ,fkCentralityMethod(0)
 127 ,fCentrality(0)
 128 ,fCentralityMin(0)
 129 ,fCentralityMax(0)
 130 ,fInvmassCut(0)
 131 ,fpTCut(0)
 132 ,fTrigger(0)
 133 ,fPhi(0)
 134 ,fEta(0)
 135 ,fVZEROA(0)
 136 ,fVZEROC(0)
 137 ,fTPCM(0)
 138 ,fNoEvents(0)
 139 ,fTrkEovPBef(0)
 140 //,fdEdxBef(0)
 141 ,fInclusiveElecPt(0)
 142 ,fTPCnsigma(0)
 143 ,fTPCnsigmaAft(0)
 144 ,fCentralityPass(0)
 145 ,fCentralityNoPass(0)
 146 ,fInvmassLS1(0)
 147 ,fInvmassULS1(0)
 148 ,fPhotoElecPt(0)
 149 ,fSemiInclElecPt(0)
 150 ,fULSElecPt(0)
 151 ,fLSElecPt(0)
 152 ,fminTPC(-1)
 153 ,fmaxTPC(3)
 154 ,fminEovP(0.8)
 155 ,fmaxEovP(1.2)
 156 ,fminM20(0.03)
 157 ,fmaxM20(0.3)
 158 ,fminM02(0.03)
 159 ,fmaxM02(0.5)
 160 ,fDispersion(1)
 161 ,fMultCorAfterCuts(0)
 162 ,fMultvsCentr(0)
 163 ,fSubEventDPhiv2(0)
 164 ,EPVzA(0)
 165 ,EPVzC(0)
 166 ,EPTPC(0)
 167 ,fV2Phi(0)
 168 ,fSparseElectronHadron(0)
 169 ,fvertex(0)
 170 ,fMultCorBeforeCuts(0)
 171 ,fSideBandsFlow(kFALSE)
 172 ,fPhiminusPsi(kFALSE)
 173 ,fFlowEventCont(0) //! flow events (one for each inv mass band)
 174 ,fpurity(kFALSE)
 175 ,fSparseElectronpurity(0)
 176 ,fOpeningAngleLS(0)
 177 ,fOpeningAngleULS(0)
 178 ,fNonHFE(new AliSelectNonHFE)
 179 ,fDCA(0)
 180 ,fOpeningAngleCut(0)
 181 ,fOP_angle(0)
 182 ,fAssoTPCCluster(0)
 183 ,fAssoITSRefit(0)
 184 ,fMultCut(0)
 185 ,fMultCorAfterCentrBeforeCuts(0)
 186 ,fMultCorAfterVZTRKComp(0)
 187 ,fCentralityBeforePileup(0)
 188 ,fCentralityAfterVZTRK(0)
 189 ,fCentralityAfterCorrCut(0)
 190 ,fMultCorAfterCorrCut(0)
 191 ,EPVz(0)
 192 ,EPTPCp(0)
 193 ,EPTPCn(0)
 194 ,fSubEventDPhiv2new(0)
 195 ,fV2Phivzerotot(0)
 196 ,fHistCentrDistr(0x0)
 197 ,fCentralityNoPassForFlattening(0)
 198 ,fInvmassLS1highpt(0)
 199 ,fInvmassULS1highpt(0)
 200 ,fSparsephipsiULS(0)
 201 ,fSparsephipsiLS(0)
 202 ,fSparseMassULS(0)
 203 ,fSparseMassLS(0)
 204 {
 205     //Named constructor
 206
 207     fPID = new AliHFEpid("hfePid");
 208     // Define input and output slots here
 209     // Input slot #0 works with a TChain
 210     DefineInput(0, TChain::Class());
 211     // Output slot #0 id reserved by the base class for AOD
 212     // Output slot #1 writes into a TH1 container
 213     // DefineOutput(1, TH1I::Class());
 214     DefineOutput(1, TList::Class());
 215     DefineOutput(2, AliFlowEventSimple::Class());
 216     if(fSideBandsFlow){
 217         DefineOutput(3, AliFlowEventSimple::Class());
 218     }
 219     //  DefineOutput(3, TTree::Class());
 220 }
 221
 222 //________________________________________________________________________
 223 AliAnalysisTaskFlowTPCEMCalQCSP::AliAnalysisTaskFlowTPCEMCalQCSP()
 224 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskFlowTPCEMCalQCSP")
 225 ,fDebug(0)
 226 ,fAOD(0)
 227 ,fGeom(0)
 228 ,fOutputList(0)
 229 ,fCuts(0)
 230 ,fIdentifiedAsOutInz(kFALSE)
 231 ,fPassTheEventCut(kFALSE)
 232 ,fCFM(0)
 233 ,fPID(0)
 234 ,fPIDqa(0)
 235 ,fCutsRP(0)     // track cuts for reference particles
 236 ,fNullCuts(0) // dummy cuts for flow event tracks
 237 ,fFlowEvent(0) //! flow events (one for each inv mass band)
 238 ,fkCentralityMethod(0)
 239 ,fCentrality(0)
 240 ,fCentralityMin(0)
 241 ,fCentralityMax(0)
 242 ,fInvmassCut(0)
 243 ,fpTCut(0)
 244 ,fTrigger(0)
 245 ,fPhi(0)
 246 ,fEta(0)
 247 ,fVZEROA(0)
 248 ,fVZEROC(0)
 249 ,fTPCM(0)
 250 ,fNoEvents(0)
 251 ,fTrkEovPBef(0)
 252 //,fdEdxBef(0)
 253 ,fInclusiveElecPt(0)
 254 ,fTPCnsigma(0)
 255 ,fTPCnsigmaAft(0)
 256 ,fCentralityPass(0)
 257 ,fCentralityNoPass(0)
 258 ,fInvmassLS1(0)
 259 ,fInvmassULS1(0)
 260 ,fPhotoElecPt(0)
 261 ,fSemiInclElecPt(0)
 262 ,fULSElecPt(0)
 263 ,fLSElecPt(0)
 264 ,fminTPC(-1)
 265 ,fmaxTPC(3)
 266 ,fminEovP(0.8)
 267 ,fmaxEovP(1.2)
 268 ,fminM20(0.03)
 269 ,fmaxM20(0.3)
 270 ,fminM02(0.03)
 271 ,fmaxM02(0.5)
 272 ,fDispersion(1)
 273 ,fMultCorAfterCuts(0)
 274 ,fMultvsCentr(0)
 275 ,fSubEventDPhiv2(0)
 276 ,EPVzA(0)
 277 ,EPVzC(0)
 278 ,EPTPC(0)
 279 ,fV2Phi(0)
 280 ,fSparseElectronHadron(0)
 281 ,fvertex(0)
 282 ,fMultCorBeforeCuts(0)
 283 ,fSideBandsFlow(kFALSE)
 284 ,fPhiminusPsi(kFALSE)
 285 ,fFlowEventCont(0) //! flow events (one for each inv mass band)
 286 ,fpurity(kFALSE)
 287 ,fSparseElectronpurity(0)
 288 ,fOpeningAngleLS(0)
 289 ,fOpeningAngleULS(0)
 290 ,fNonHFE(new AliSelectNonHFE)
 291 ,fDCA(0)
 292 ,fOpeningAngleCut(0)
 293 ,fOP_angle(0)
 294 ,fAssoTPCCluster(0)
 295 ,fAssoITSRefit(0)
 296 ,fMultCut(0)
 297 ,fMultCorAfterCentrBeforeCuts(0)
 298 ,fMultCorAfterVZTRKComp(0)
 299 ,fCentralityBeforePileup(0)
 300 ,fCentralityAfterVZTRK(0)
 301 ,fCentralityAfterCorrCut(0)
 302 ,fMultCorAfterCorrCut(0)
 303 ,EPVz(0)
 304 ,EPTPCp(0)
 305 ,EPTPCn(0)
 306 ,fSubEventDPhiv2new(0)
 307 ,fV2Phivzerotot(0)
 308 ,fHistCentrDistr(0x0)
 309 ,fCentralityNoPassForFlattening(0)
 310 ,fInvmassLS1highpt(0)
 311 ,fInvmassULS1highpt(0)
 312 ,fSparsephipsiULS(0)
 313 ,fSparsephipsiLS(0)
 314 ,fSparseMassULS(0)
 315 ,fSparseMassLS(0)
 316 {
 317     //Default constructor
 318     fPID = new AliHFEpid("hfePid");
 319     // Constructor
 320     // Define input and output slots here
 321     // Input slot #0 works with a TChain
 322     DefineInput(0, TChain::Class());
 323     // Output slot #0 id reserved by the base class for AOD
 324     // Output slot #1 writes into a TH1 container
 325     // DefineOutput(1, TH1I::Class());
 326     DefineOutput(1, TList::Class());
 327     DefineOutput(2, AliFlowEventSimple::Class());
 328     //  DefineOutput(3, TTree::Class());
 329     if(fSideBandsFlow){
 330         DefineOutput(3, AliFlowEventSimple::Class());
 331     }
 332     //DefineOutput(3, TTree::Class());
 333 }
 334 //_________________________________________
 335
 336 AliAnalysisTaskFlowTPCEMCalQCSP::~AliAnalysisTaskFlowTPCEMCalQCSP()
 337 {
 338     //Destructor
 339
 340     delete fOutputList;
 341     delete fGeom;
 342     delete fPID;
 343     delete fCFM;
 344     delete fPIDqa;
 345     if (fDCA)  delete fNonHFE;
 346     if (fOutputList) delete fOutputList;
 347     if (fFlowEvent) delete fFlowEvent;
 348     if (fFlowEventCont) delete fFlowEventCont;
 349
 350 }
 351 //_________________________________________
 352
 353 void AliAnalysisTaskFlowTPCEMCalQCSP::UserExec(Option_t*)
 354 {
 355     //Main loop
 356     //Called for each event
 357
 358     // create pointer to event
 359
 360     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 361     if (!fAOD)
 362     {
 363         printf("ERROR: fAOD not available\n");
 364         return;
 365     }
 366
 367     if(!fCuts)
 368     {
 369         AliError("HFE cuts not available");
 370         return;
 371     }
 372
 373     if(!fPID->IsInitialized())
 374     {
 375         // Initialize PID with the given run number
 376         AliWarning("PID not initialised, get from Run no");
 377         fPID->InitializePID(fAOD->GetRunNumber());
 378     }
 379
 380     //  cout << "kTrigger   ==   " << fTrigger <<endl;
 381
 382     if(fTrigger==0){
 383         if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral) ) return;
 384     }
 385     if(fTrigger==1){
 386         if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
 387     }
 388     if(fTrigger==2){
 389         if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA) ) return;
 390     }
 391     if(fTrigger==3){
 392         if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) ) return;
 393     }
 394     if(fTrigger==4){
 395         if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kCentral | AliVEvent::kSemiCentral))) return;
 396     }
 397
 398
 399     //---------------CENTRALITY AND EVENT SELECTION-----------------------
 400     Int_t fNOtrks =  fAOD->GetNumberOfTracks();
 401     Float_t vtxz = -999;
 402     const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
 403     if (!trkVtx || trkVtx->GetNContributors()<=0)return;
 404     TString vtxTtl = trkVtx->GetTitle();
 405     if (!vtxTtl.Contains("VertexerTracks"))return;
 406     const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
 407     if (!spdVtx || spdVtx->GetNContributors()<=0)return;
 408     if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
 409     vtxz = trkVtx->GetZ();
 410     if(TMath::Abs(vtxz)>10)return;
 411
 412     // Event cut
 413     if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
 414     if(fNOtrks<2) return;
 415
 416     Bool_t pass = kFALSE; //to select centrality and pile up protection
 417     CheckCentrality(fAOD,pass);
 418     if(!pass)return;
 419     fvertex->Fill(vtxz);
 420
 421     fNoEvents->Fill(0);
 422     PlotVZeroMultiplcities(fAOD);
 423
 424     SetNullCuts(fAOD);
 425     PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent);    //Calculate event plane Qvector and EP resolution for inclusive
 426
 427     if(fSideBandsFlow){
 428         PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEventCont);    //Calculate event plane Qvector and EP resolution for inclusive
 429     }
 430
 431     AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 432     if(!pidResponse)
 433     {
 434         AliDebug(1, "Using default PID Response");
 435         pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 436     }
 437
 438     fPID->SetPIDResponse(pidResponse);
 439     fCFM->SetRecEventInfo(fAOD);
 440
 441     // Look for kink mother
 442     Int_t numberofvertices = fAOD->GetNumberOfVertices();
 443     Double_t listofmotherkink[numberofvertices];
 444     Int_t numberofmotherkink = 0;
 445     for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
 446         AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
 447         if(!aodvertex) continue;
 448         if(aodvertex->GetType()==AliAODVertex::kKink) {
 449             AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
 450             if(!mother) continue;
 451             Int_t idmother = mother->GetID();
 452             listofmotherkink[numberofmotherkink] = idmother;
 453             //printf("ID %d\n",idmother);
 454             numberofmotherkink++;
 455         }
 456     }
 457
 458
 459     //=============================================V0EP from Alex======================================================================
 460     Double_t qxEPa = 0, qyEPa = 0;
 461     Double_t qxEPc = 0, qyEPc = 0;
 462     Double_t qxEP = 0, qyEP = 0;
 463
 464     Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
 465     Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
 466     Double_t evPlAngV0 = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 10, 2, qxEP, qyEP);
 467
 468
 469     Double_t Qx2 = 0, Qy2 = 0;
 470     Double_t Qx2p = 0, Qy2p = 0;
 471     Double_t Qx2n = 0, Qy2n = 0;
 472
 473     for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
 474
 475         AliAODTrack* aodTrack = fAOD->GetTrack(iT);
 476
 477         if (!aodTrack)
 478             continue;
 479
 480         if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
 481             continue;
 482
 483         if (!aodTrack->TestFilterBit(128))
 484             continue;
 485
 486
 487         if(aodTrack->Eta()>0 && aodTrack->Eta()<0.8){
 488
 489             Qx2p += TMath::Cos(2*aodTrack->Phi());
 490             Qy2p += TMath::Sin(2*aodTrack->Phi());
 491         }
 492         if(aodTrack->Eta()<0 && aodTrack->Eta()> -0.8){
 493
 494             Qx2n += TMath::Cos(2*aodTrack->Phi());
 495             Qy2n += TMath::Sin(2*aodTrack->Phi());
 496         }
 497
 498
 499         Qx2 += TMath::Cos(2*aodTrack->Phi());
 500         Qy2 += TMath::Sin(2*aodTrack->Phi());
 501
 502
 503
 504
 505     }
 506
 507     Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
 508     Double_t evPlAngTPCn = TMath::ATan2(Qy2n, Qx2n)/2.;
 509     Double_t evPlAngTPCp = TMath::ATan2(Qy2p, Qx2p)/2.;
 510
 511     EPVzA->Fill(evPlAngV0A);
 512     EPVzC->Fill(evPlAngV0C);
 513     EPTPC->Fill(evPlAngTPC);
 514
 515     EPTPCn->Fill(evPlAngTPCn);
 516     EPTPCp->Fill(evPlAngTPCp);
 517     EPVz->Fill(evPlAngV0);
 518
 519
 520
 521
 522     fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
 523     fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
 524     fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
 525
 526
 527     fSubEventDPhiv2new->Fill(0.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCp))); // vzero - tpcp
 528     fSubEventDPhiv2new->Fill(1.5, TMath::Cos(2.*(evPlAngV0-evPlAngTPCn))); // vzero - tpcn
 529     fSubEventDPhiv2new->Fill(2.5, TMath::Cos(2.*(evPlAngTPCp-evPlAngTPCn))); // tpcp - tpcn
 530
 531
 532
 533
 534     //====================================================================================================================
 535
 536
 537     AliAODTrack *track = NULL;
 538
 539     // Track loop
 540     for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
 541     {
 542         track = fAOD->GetTrack(iTracks);
 543         if (!track)
 544         {
 545             printf("ERROR: Could not receive track %d\n", iTracks);
 546             continue;
 547         }
 548
 549         if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
 550         //----------hfe begin---------
 551         if(track->Eta()<-0.7 || track->Eta()>0.7)    continue;    //eta cuts on candidates
 552
 553         // RecKine: ITSTPC cuts
 554         if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 555
 556         // Reject kink mother
 557         Bool_t kinkmotherpass = kTRUE;
 558         for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
 559             if(track->GetID() == listofmotherkink[kinkmother]) {
 560                 kinkmotherpass = kFALSE;
 561                 continue;
 562             }
 563         }
 564         if(!kinkmotherpass) continue;
 565
 566         // RecPrim
 567         //  if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;  //deleted for DCA absence
 568         // HFEcuts: ITS layers cuts
 569         if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 570         // HFE cuts: TPC PID cleanup
 571         if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 572
 573         Double_t fClsE = -999, p = -999, fEovP=-999, pt = -999, fTPCnSigma=0;
 574         // Track extrapolation
 575         Int_t fClsId = track->GetEMCALcluster();
 576         if(fClsId < 0) continue;
 577         AliAODCaloCluster *cluster = fAOD->GetCaloCluster(fClsId);
 578         if(TMath::Abs(cluster->GetTrackDx()) > 0.05 || TMath::Abs(cluster->GetTrackDz()) > 0.05) continue;
 579
 580         pt = track->Pt();         //pt track after cuts
 581         if(pt<fpTCut) continue;
 582         fClsE = cluster->E();
 583         p = track->P();
 584         // dEdx = track->GetTPCsignal();
 585         fEovP = fClsE/p;
 586         fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000;
 587         Double_t m20 =cluster->GetM20();
 588         Double_t m02 =cluster->GetM02();
 589         Double_t disp=cluster->GetDispersion();
 590         if(fTPCnSigma >= -1 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,fEovP);
 591         fTPCnsigma->Fill(p,fTPCnSigma);
 592         //   fdEdxBef->Fill(p,dEdx);
 593         Double_t eta = track->Eta();
 594         Double_t phi = track->Phi();
 595         //-----------------------Phiminupsi method to remove the contamination-----------------------------------------------
 596         //-----------------------fTPCnSigma < -3.5 hadrons will be selected from this region--------------------------
 597         Float_t dPhi_aeh = TVector2::Phi_0_2pi(phi - evPlAngV0A);
 598         if(dPhi_aeh > TMath::Pi()) dPhi_aeh = dPhi_aeh - TMath::Pi();
 599         Float_t dPhi_ceh = TVector2::Phi_0_2pi(phi - evPlAngV0C);
 600         if(dPhi_ceh > TMath::Pi()) dPhi_ceh = dPhi_ceh - TMath::Pi();
 601
 602         if(fPhiminusPsi){
 603             Double_t valueElh[8] = {
 604                 pt,
 605                 fEovP,
 606                 fTPCnSigma,
 607                 m20,
 608                 m02,
 609                 disp,
 610                 dPhi_aeh,
 611                 dPhi_ceh};
 612             fSparseElectronHadron->Fill(valueElh);
 613         }
 614         //----------------------------------------------------------------------------------------------------------
 615         //---------------------------From here usual electron selection---------------------------------------------
 616         //----------------------------------------------------------------------------------------------------------
 617         if(m20 < fminM20 || m20 > fmaxM20) continue;
 618         if(m02 < fminM02 || m02 > fmaxM02) continue;
 619         if(disp > fDispersion ) continue;
 620         //---------------------------------for purity---------------------------------------------------------------
 621         if(fpurity){
 622             Double_t valuepurity[3] = {
 623                 pt,
 624                 fEovP,
 625                 fTPCnSigma};
 626             fSparseElectronpurity->Fill(valuepurity);
 627         }
 628         //----------------------------------------------------------------------------------------------------------
 629         //----------------------------------------------------------------------------------------------------------
 630         if(fTPCnSigma < fminTPC || fTPCnSigma > fmaxTPC) continue; //cuts on nsigma tpc and EoP
 631         //===============================Flow Event for Contamination=============================================
 632         if(fSideBandsFlow){
 633             if(fEovP>0 && fEovP<0.6){
 634                 AliFlowTrack *sTrackCont = new AliFlowTrack();
 635                 sTrackCont->Set(track);
 636                 sTrackCont->SetID(track->GetID());
 637                 sTrackCont->SetForRPSelection(kTRUE);
 638                 sTrackCont->SetForPOISelection(kTRUE);
 639                 sTrackCont->SetMass(2637);
 640                 for(int iRPs=0; iRPs!=fFlowEventCont->NumberOfTracks(); ++iRPs)
 641                 {
 642                     //   cout << " no of rps " << iRPs << endl;
 643                     AliFlowTrack *iRPCont = dynamic_cast<AliFlowTrack*>(fFlowEventCont->GetTrack( iRPs ));
 644                     if (!iRPCont) continue;
 645                     if (!iRPCont->InRPSelection()) continue;
 646                     if( sTrackCont->GetID() == iRPCont->GetID())
 647                     {
 648                         if(fDebug) printf(" was in RP set");
 649                         //       cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set" <<endl;
 650                         iRPCont->SetForRPSelection(kFALSE);
 651                       //  fFlowEventCont->SetNumberOfRPs(fFlowEventCont->GetNumberOfRPs() - 1);
 652                     }
 653                 } //end of for loop on RPs
 654                 fFlowEventCont->InsertTrack(((AliFlowTrack*) sTrackCont));
 655                 fFlowEventCont->SetNumberOfPOIs(fFlowEventCont->GetNumberOfPOIs()+1);
 656
 657             }
 658         }
 659         //==========================================================================================================
 660         //===============================From here eovP cut is used fro QC, SP and EPV0=============================
 661         if(fEovP < fminEovP || fEovP >fmaxEovP) continue;
 662         //==========================================================================================================
 663         //============================Event Plane Method with V0====================================================
 664         Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
 665         Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
 666         Double_t v2Phi[3] = {
 667             v2PhiV0A,
 668             v2PhiV0C,
 669             pt};
 670         fV2Phi->Fill(v2Phi);
 671
 672         Double_t v2PhiVz = TMath::Cos(2*(phi - evPlAngV0));
 673         Double_t v2PhiV0tot[2] = {
 674             v2PhiVz,
 675             pt};
 676         fV2Phivzerotot->Fill(v2PhiV0tot);
 677         //=========================================================================================================
 678         fTPCnsigmaAft->Fill(p,fTPCnSigma);
 679         fInclusiveElecPt->Fill(pt);
 680         fPhi->Fill(phi);
 681         fEta->Fill(eta);
 682         //----------------------Flow of Inclusive Electrons--------------------------------------------------------
 683         AliFlowTrack *sTrack = new AliFlowTrack();
 684         sTrack->Set(track);
 685         sTrack->SetID(track->GetID());
 686         sTrack->SetForRPSelection(kTRUE);
 687         sTrack->SetForPOISelection(kTRUE);
 688         sTrack->SetMass(263732);
 689         for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
 690         {
 691             //   cout << " no of rps " << iRPs << endl;
 692             AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
 693             if (!iRP) continue;
 694             if (!iRP->InRPSelection()) continue;
 695             if( sTrack->GetID() == iRP->GetID())
 696             {
 697                 if(fDebug) printf(" was in RP set");
 698                 //       cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set" <<endl;
 699                 iRP->SetForRPSelection(kFALSE);
 700                // fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
 701             }
 702         } //end of for loop on RPs
 703         fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
 704         fFlowEvent->SetNumberOfPOIs(fFlowEvent->GetNumberOfPOIs()+1);
 705
 706
 707         if(fDCA){
 708             //----------------------Selection of Photonic Electrons DCA-----------------------------
 709             fNonHFE = new AliSelectNonHFE();
 710             fNonHFE->SetAODanalysis(kTRUE);
 711             fNonHFE->SetInvariantMassCut(fInvmassCut);
 712             if(fOP_angle) fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
 713             //fNonHFE->SetChi2OverNDFCut(fChi2Cut);
 714             //if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
 715             fNonHFE->SetAlgorithm("DCA"); //KF
 716             fNonHFE->SetPIDresponse(pidResponse);
 717             fNonHFE->SetTrackCuts(-3,3);
 718
 719             fNonHFE->SetHistAngleBack(fOpeningAngleLS);
 720             fNonHFE->SetHistAngle(fOpeningAngleULS);
 721             //fNonHFE->SetHistDCABack(fDCABack);
 722             //fNonHFE->SetHistDCA(fDCA);
 723             fNonHFE->SetHistMassBack(fInvmassLS1);
 724             fNonHFE->SetHistMass(fInvmassULS1);
 725
 726             fNonHFE->FindNonHFE(iTracks,track,fAOD);
 727
 728             // Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
 729             // Int_t *fLsPartner = fNonHFE->GetPartnersLS();
 730             // Bool_t fUlsIsPartner = kFALSE;
 731             // Bool_t fLsIsPartner = kFALSE;
 732             if(fNonHFE->IsULS()){
 733                 for(Int_t kULS =0; kULS < fNonHFE->GetNULS(); kULS++){
 734                     fULSElecPt->Fill(track->Pt());
 735                 }
 736             }
 737
 738             if(fNonHFE->IsLS()){
 739                 for(Int_t kLS =0; kLS < fNonHFE->GetNLS(); kLS++){
 740                     fLSElecPt->Fill(track->Pt());
 741                 }
 742             }
 743         }
 744
 745         if(!fDCA){
 746             //----------------------Selection of Photonic Electrons KFParticle-----------------------------
 747             Bool_t fFlagPhotonicElec = kFALSE;
 748             SelectPhotonicElectron(iTracks,track,fEovP, evPlAngV0, fFlagPhotonicElec);
 749             if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
 750             // Semi inclusive electron
 751             if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
 752         }
 753
 754
 755
 756     }//end loop on track
 757
 758     PostData(1, fOutputList);
 759     PostData(2, fFlowEvent);
 760     if(fSideBandsFlow){
 761         PostData(3, fFlowEventCont);
 762     }
 763
 764     //----------hfe end---------
 765 }
 766 //_________________________________________
 767 void AliAnalysisTaskFlowTPCEMCalQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track,Double_t fEovP,Double_t evPlAngV0, Bool_t &fFlagPhotonicElec)
 768 {
 769     //Identify non-heavy flavour electrons using Invariant mass method KF
 770
 771     Bool_t flagPhotonicElec = kFALSE;
 772
 773     for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
 774         AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
 775         if (!trackAsso) {
 776             printf("ERROR: Could not receive track %d\n", jTracks);
 777             continue;
 778         }
 779         //  if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
 780         if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
 781         //    if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit) || (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
 782
 783         if(fAssoITSRefit){
 784             if(!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)) continue;
 785         }
 786
 787         if(!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)) continue;
 788
 789         if(jTracks == itrack) continue;
 790         Double_t ptAsso=-999., nsigma=-999.0;
 791         Double_t mass=-999., width = -999;
 792         Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 793         Double_t openingAngle = -999.;
 794         Double_t ptcutonmasshighpt = track->Pt();
 795
 796         ptAsso = trackAsso->Pt();
 797         Short_t chargeAsso = trackAsso->Charge();
 798         Short_t charge = track->Charge();
 799
 800         nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
 801
 802         //80
 803         if(trackAsso->GetTPCNcls() < fAssoTPCCluster) continue;
 804         if(nsigma < -3 || nsigma > 3) continue;
 805         if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
 806         if(ptAsso <0.3) continue;
 807
 808         Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 809         if(charge>0) fPDGe1 = -11;
 810         if(chargeAsso>0) fPDGe2 = -11;
 811
 812         if(charge == chargeAsso) fFlagLS = kTRUE;
 813         if(charge != chargeAsso) fFlagULS = kTRUE;
 814
 815         AliKFParticle::SetField(fAOD->GetMagneticField());
 816         AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
 817         AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
 818         AliKFParticle recg(ge1, ge2);
 819
 820         if(recg.GetNDF()<1) continue;
 821         Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 822         if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 823         recg.GetMass(mass,width);
 824
 825         openingAngle = ge1.GetAngle(ge2);
 826         if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
 827         if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
 828         if(fOP_angle)if(openingAngle > fOpeningAngleCut) continue;
 829
 830
 831         if(fFlagLS) fInvmassLS1->Fill(mass);
 832         if(fFlagULS) fInvmassULS1->Fill(mass);
 833
 834         if(fFlagULS){
 835         Double_t MassSparseULS[3] = {
 836                 track->Pt(),
 837                 mass
 838                 };
 839                 fSparseMassULS->Fill(MassSparseULS);
 840          }
 841         if(fFlagLS){
 842         Double_t MassSparseLS[3] = {
 843                 track->Pt(),
 844                 mass
 845                 };
 846                 fSparseMassLS->Fill(MassSparseLS);
 847          }
 848
 849
 850         if(ptcutonmasshighpt >= 8.){
 851             if(fFlagLS) fInvmassLS1highpt->Fill(mass);
 852             if(fFlagULS) fInvmassULS1highpt->Fill(mass);
 853         }
 854
 855
 856         if(mass<fInvmassCut){
 857             if(fFlagULS){fULSElecPt->Fill(track->Pt());}
 858             if(fFlagLS){fLSElecPt->Fill(track->Pt());}
 859         }
 860
 861
 862
 863         Double_t phi = track->Phi();
 864         Float_t DeltaPhi_eEP = TVector2::Phi_0_2pi(phi - evPlAngV0);
 865         if(DeltaPhi_eEP > TMath::Pi()) {DeltaPhi_eEP = DeltaPhi_eEP - TMath::Pi();}
 866
 867
 868  if(mass<fInvmassCut){
 869         if(fFlagULS){
 870                 Double_t ulsSparse[3] = {
 871                 track->Pt(),
 872                 fEovP,
 873                 DeltaPhi_eEP
 874                 };
 875                 fSparsephipsiULS->Fill(ulsSparse);
 876                 }
 877         if(fFlagLS){
 878                 Double_t lsSparse[3] = {
 879                 track->Pt(),
 880                 fEovP,
 881                 DeltaPhi_eEP
 882                 };
 883                 fSparsephipsiLS->Fill(lsSparse);
 884                 }
 885         }
 886
 887
 888
 889         if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
 890             flagPhotonicElec = kTRUE;
 891         }
 892     }//track loop
 893     fFlagPhotonicElec = flagPhotonicElec;
 894 }
 895 //__________________________________________________________________________________
 896 void AliAnalysisTaskFlowTPCEMCalQCSP::UserCreateOutputObjects()
 897 {
 898     //Create histograms
 899
 900     //----------hfe initialising begin---------
 901     fNullCuts = new AliFlowTrackCuts("null_cuts");
 902
 903     AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 904     cc->SetNbinsMult(10000);
 905     cc->SetMultMin(0);
 906     cc->SetMultMax(10000);
 907
 908     cc->SetNbinsPt(100);
 909     cc->SetPtMin(0);
 910     cc->SetPtMax(50);
 911
 912     cc->SetNbinsPhi(180);
 913     cc->SetPhiMin(0.0);
 914     cc->SetPhiMax(TMath::TwoPi());
 915
 916     cc->SetNbinsEta(30);
 917     cc->SetEtaMin(-7.0);
 918     cc->SetEtaMax(+7.0);
 919
 920     cc->SetNbinsQ(500);
 921     cc->SetQMin(0.0);
 922     cc->SetQMax(3.0);
 923
 924     //--------Initialize PID
 925     fPID->SetHasMCData(kFALSE);
 926     if(!fPID->GetNumberOfPIDdetectors())
 927     {
 928         fPID->AddDetector("TPC", 0);
 929         fPID->AddDetector("EMCAL", 1);
 930     }
 931
 932     fPID->SortDetectors();
 933     fPIDqa = new AliHFEpidQAmanager();
 934     fPIDqa->Initialize(fPID);
 935
 936     //--------Initialize correction Framework and Cuts
 937     fCFM = new AliCFManager;
 938     const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 939     fCFM->SetNStepParticle(kNcutSteps);
 940     for(Int_t istep = 0; istep < kNcutSteps; istep++)
 941         fCFM->SetParticleCutsList(istep, NULL);
 942
 943     if(!fCuts){
 944         AliWarning("Cuts not available. Default cuts will be used");
 945         fCuts = new AliHFEcuts;
 946         fCuts->CreateStandardCuts();
 947     }
 948
 949     fCuts->SetAOD();
 950     fCuts->Initialize(fCFM);
 951     //----------hfe initialising end--------
 952     //---------Output Tlist
 953     fOutputList = new TList();
 954     fOutputList->SetOwner();
 955     fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 956
 957     fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
 958     fOutputList->Add(fNoEvents);
 959
 960     fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",1000,0,50,200,-10,10);
 961     fOutputList->Add(fTPCnsigma);
 962
 963     fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",1000,0,50,200,-10,10);
 964     fOutputList->Add(fTPCnsigmaAft);
 965
 966     fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",1000,0,50,100,0,2);
 967     fOutputList->Add(fTrkEovPBef);
 968
 969     //    fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",1000,0,50,150,0,150);
 970     //     fOutputList->Add(fdEdxBef);
 971
 972     fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",1000,0,100);
 973     fOutputList->Add(fInclusiveElecPt);
 974
 975     fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",1000,0,100);
 976     fOutputList->Add(fPhotoElecPt);
 977
 978     fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",1000,0,100);
 979     fOutputList->Add(fSemiInclElecPt);
 980
 981     fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",1000,0,100);
 982     fOutputList->Add(fULSElecPt);
 983
 984     fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",1000,0,100);
 985     fOutputList->Add(fLSElecPt);
 986
 987     fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
 988     fOutputList->Add(fInvmassLS1);
 989
 990     fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
 991     fOutputList->Add(fInvmassULS1);
 992
 993     fInvmassLS1highpt = new TH1F("fInvmassLS1highpt", "Inv mass of LS (e,e); mass(GeV/c^2) highpt; counts;", 1000,0,1.0);
 994     fOutputList->Add(fInvmassLS1highpt);
 995
 996     fInvmassULS1highpt = new TH1F("fInvmassULS1highpt", "Inv mass of ULS (e,e); mass(GeV/c^2) highpt; counts;", 1000,0,1.0);
 997     fOutputList->Add(fInvmassULS1highpt);
 998
 999     fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
1000     fOutputList->Add(fCentralityPass);
1001
1002     fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
1003     fOutputList->Add(fCentralityNoPass);
1004
1005     fCentralityNoPassForFlattening = new TH1F("fCentralityNoPassForFlattening", "Centrality No Pass for flattening", 101, -1, 100);
1006     fOutputList->Add(fCentralityNoPassForFlattening);
1007
1008     fCentralityBeforePileup = new TH1F("fCentralityBeforePileup", "fCentralityBeforePileup Pass", 101, -1, 100);
1009     fOutputList->Add(fCentralityBeforePileup);
1010
1011     fCentralityAfterVZTRK = new TH1F("fCentralityAfterVZTRK", "fCentralityAfterVZTRK Pass", 101, -1, 100);
1012     fOutputList->Add(fCentralityAfterVZTRK);
1013
1014     fCentralityAfterCorrCut = new TH1F("fCentralityAfterCorrCut", "fCentralityAfterCorrCut Pass", 101, -1, 100);
1015     fOutputList->Add(fCentralityAfterCorrCut);
1016
1017     fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
1018     fOutputList->Add(fPhi);
1019
1020     fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
1021     fOutputList->Add(fEta);
1022
1023     fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
1024     fOutputList->Add(fVZEROA);
1025
1026     fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
1027     fOutputList->Add(fVZEROC);
1028
1029     fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
1030     fOutputList->Add(fTPCM);
1031
1032     fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
1033     fOutputList->Add(fvertex);
1034
1035     fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1036     fOutputList->Add(fMultCorBeforeCuts);
1037
1038     fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1039     fOutputList->Add(fMultCorAfterCuts);
1040
1041     fMultCorAfterCentrBeforeCuts = new TH2F("fMultCorAfterCentrBeforeCuts", "TPC vs Global multiplicity (After CC before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1042     fOutputList->Add(fMultCorAfterCentrBeforeCuts);
1043
1044     fMultCorAfterVZTRKComp = new TH2F("fMultCorAfterVZTRKComp", "TPC vs Global multiplicity (After V0-TRK); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1045     fOutputList->Add(fMultCorAfterVZTRKComp);
1046
1047     fMultCorAfterCorrCut = new TH2F("fMultCorAfterCorrCut", "TPC vs Global multiplicity (After CorrCut); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
1048     fOutputList->Add(fMultCorAfterCorrCut);
1049
1050     fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
1051     fOutputList->Add(fMultvsCentr);
1052
1053     fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1054     fOutputList->Add(fOpeningAngleLS);
1055
1056     fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1057     fOutputList->Add(fOpeningAngleULS);
1058
1059
1060     //----------------------------------------------------------------------------
1061     EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
1062     fOutputList->Add(EPVzA);
1063     EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
1064     fOutputList->Add(EPVzC);
1065     EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
1066     fOutputList->Add(EPTPC);
1067
1068
1069     EPVz = new TH1D("EPVz", "EPVz", 80, -2, 2);
1070     fOutputList->Add(EPVz);
1071     EPTPCp = new TH1D("EPTPCp", "EPTPCp", 80, -2, 2);
1072     fOutputList->Add(EPTPCp);
1073     EPTPCn = new TH1D("EPTPCn", "EPTPCn", 80, -2, 2);
1074     fOutputList->Add(EPTPCn);
1075
1076     //----------------------------------------------------------------------------
1077     fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
1078     fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1079     fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1080     fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1081     fOutputList->Add(fSubEventDPhiv2);
1082
1083
1084
1085     fSubEventDPhiv2new = new TProfile("fSubEventDPhiv2new", "fSubEventDPhiv2new", 3, 0, 3);
1086     fSubEventDPhiv2new->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1087     fSubEventDPhiv2new->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1088     fSubEventDPhiv2new->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1089     fOutputList->Add(fSubEventDPhiv2new);
1090     //================================Event Plane with VZERO A & C=====================
1091     const Int_t nPtBins = 12;
1092     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, 10, 13};
1093     // v2A, v2C, pt
1094     Int_t    bins[3] = {  50,  50, nPtBins};
1095     Double_t xmin[3] = { -1., -1.,   0};
1096     Double_t xmax[3] = {  1.,  1.,   13};
1097     fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
1098     // Set bin limits for axes which are not standard binned
1099     fV2Phi->SetBinEdges(2, binsPt);
1100     // set axes titles
1101     fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
1102     fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
1103     fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1104     fOutputList->Add(fV2Phi);
1105
1106
1107
1108     //================================Event Plane with VZERO=====================
1109     // const Int_t nPtBins = 10;
1110     // 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};
1111     // v2, pt
1112     Int_t    binsV[2] = {  50,  nPtBins};
1113     Double_t xminV[2] = { -1.,   0};
1114     Double_t xmaxV[2] = {  1.,   13};
1115     fV2Phivzerotot = new THnSparseF("fV2Phivzerotot", "v2:pt", 2, binsV, xminV, xmaxV);
1116     // Set bin limits for axes which are not standard binned
1117     fV2Phivzerotot->SetBinEdges(1, binsPt);
1118     // set axes titles
1119     fV2Phivzerotot->GetAxis(0)->SetTitle("v_{2} (V0)");
1120     fV2Phivzerotot->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
1121     fOutputList->Add(fV2Phivzerotot);
1122
1123
1124
1125     //----------------------------------------------------------------------------
1126
1127     if(fPhiminusPsi){
1128         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)
1129         Double_t xminvElectH[8]={0,    0, -10 ,  0,    0,    0,    0,           0};
1130         Double_t xmaxvElectH[8]={20,   2,  10 ,  2,    2,    2,  TMath::Pi(), TMath::Pi()};
1131         fSparseElectronHadron = new THnSparseD("ElectronHadron","ElectronHadron",8,binsvElectH,xminvElectH,xmaxvElectH);
1132         fSparseElectronHadron->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1133         fSparseElectronHadron->GetAxis(1)->SetTitle("EovP");
1134         fSparseElectronHadron->GetAxis(2)->SetTitle("TPCnSigma");
1135         fSparseElectronHadron->GetAxis(3)->SetTitle("M20");
1136         fSparseElectronHadron->GetAxis(4)->SetTitle("M02");
1137         fSparseElectronHadron->GetAxis(5)->SetTitle("Disp");
1138         fSparseElectronHadron->GetAxis(6)->SetTitle("phiminuspsi V0A");
1139         fSparseElectronHadron->GetAxis(7)->SetTitle("phiminuspsi V0C");
1140         fOutputList->Add(fSparseElectronHadron);
1141     }
1142     //----------------------------------------------------------------------------
1143     //----------------------------------------------------------------------------
1144     if(fpurity){
1145         Int_t binsvpurity[3]={   600,200, 200}; //pt, E/p,TPCnSigma
1146         Double_t xminvpurity[3]={0,    0, -10};
1147         Double_t xmaxvpurity[3]={30,   2,  10};
1148         fSparseElectronpurity = new THnSparseD("Electronpurity","Electronpurity",3,binsvpurity,xminvpurity,xmaxvpurity);
1149         fSparseElectronpurity->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1150         fSparseElectronpurity->GetAxis(1)->SetTitle("EovP");
1151         fSparseElectronpurity->GetAxis(2)->SetTitle("TPCnSigma");
1152         fOutputList->Add(fSparseElectronpurity);
1153     }
1154     //----------------------------------------------------------------------------
1155
1156
1157     Int_t    binsphipsi[3] = { 100,   200,           6};
1158     Double_t xminphipsi[3] = { 0.,      0,           0};
1159     Double_t xmaxphipsi[3] = { 10.,      2, TMath::Pi()};
1160     fSparsephipsiULS = new THnSparseF("fSparsephipsiULS", "pt:eop:DeltaPhiULS", 3, binsphipsi, xminphipsi, xmaxphipsi);
1161     fSparsephipsiULS->GetAxis(0)->SetTitle("pt (Gev/c)");
1162     fSparsephipsiULS->GetAxis(1)->SetTitle("eop");
1163     fSparsephipsiULS->GetAxis(2)->SetTitle("DeltaPhiULS");
1164     fOutputList->Add(fSparsephipsiULS);
1165
1166     fSparsephipsiLS = new THnSparseF("fSparsephipsiLS", "pt:eop:DeltaPhiLS", 3, binsphipsi, xminphipsi, xmaxphipsi);
1167     fSparsephipsiLS->GetAxis(0)->SetTitle("pt (Gev/c)");
1168     fSparsephipsiLS->GetAxis(1)->SetTitle("eop");
1169     fSparsephipsiLS->GetAxis(2)->SetTitle("DeltaPhiLS");
1170     fOutputList->Add(fSparsephipsiLS);
1171
1172     Int_t    binsmass[2] = { 100, 200};
1173     Double_t xminmass[2] = { 0.,  0};
1174     Double_t xmaxmass[2] = { 10., 1.};
1175     fSparseMassULS = new THnSparseF("fSparseMassULS", "pt:mass (GeV/c^{2})", 2, binsmass, xminmass, xmaxmass);
1176     fSparseMassULS->GetAxis(0)->SetTitle("pt (Gev/c)");
1177     fSparseMassULS->GetAxis(1)->SetTitle("mass");
1178     fOutputList->Add(fSparseMassULS);
1179
1180     fSparseMassLS = new THnSparseF("fSparseMassLS", "pt:mass (GeV/c^{2})", 2, binsmass, xminmass, xmaxmass);
1181     fSparseMassLS->GetAxis(0)->SetTitle("pt (Gev/c)");
1182     fSparseMassLS->GetAxis(1)->SetTitle("mass");
1183     fOutputList->Add(fSparseMassLS);
1184
1185
1186
1187     PostData(1,fOutputList);
1188     // create and post flowevent
1189     fFlowEvent = new AliFlowEvent(10000);
1190     PostData(2, fFlowEvent);
1191
1192     if(fSideBandsFlow){
1193         fFlowEventCont = new AliFlowEvent(10000);
1194         PostData(3, fFlowEventCont);
1195     }
1196 }
1197
1198 //________________________________________________________________________
1199 void AliAnalysisTaskFlowTPCEMCalQCSP::Terminate(Option_t *)
1200 {
1201     // Info("Terminate");
1202     AliAnalysisTaskSE::Terminate();
1203 }
1204 //_____________________________________________________________________________
1205 template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::PlotVZeroMultiplcities(const T* event) const
1206 {
1207     // QA multiplicity plots
1208     fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
1209     fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
1210 }
1211 //_____________________________________________________________________________
1212 template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::SetNullCuts(T* event)
1213 {
1214     //Set null cuts
1215     if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
1216     fCutsRP->SetEvent(event, MCEvent());
1217     fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
1218     fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
1219     fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
1220     fNullCuts->SetEvent(event, MCEvent());
1221 }
1222 //_____________________________________________________________________________
1223 void AliAnalysisTaskFlowTPCEMCalQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
1224 {
1225     //Prepare flow events
1226     FlowEv->ClearFast();
1227     FlowEv->Fill(fCutsRP, fNullCuts);
1228     FlowEv->SetReferenceMultiplicity(iMulti);
1229     FlowEv->DefineDeadZone(0, 0, 0, 0);
1230     //  FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
1231 }
1232 //_____________________________________________________________________________
1233 Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1234 {
1235     // Check single track cuts for a given cut step
1236     const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1237     if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1238     return kTRUE;
1239 }
1240 //_________________________________________
1241 void AliAnalysisTaskFlowTPCEMCalQCSP::CheckCentrality(AliAODEvent* event, Bool_t &centralitypass)
1242 {
1243     //============================Multiplicity TPV vs Global===============================================================================
1244     const Int_t nGoodTracks = event->GetNumberOfTracks();
1245     Float_t multTPC(0.); // tpc mult estimate
1246     Float_t multGlob(0.); // global multiplicity
1247     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
1248         AliAODTrack* trackAOD = event->GetTrack(iTracks);
1249         if (!trackAOD) continue;
1250         if (!(trackAOD->TestFilterBit(1))) continue;
1251         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;
1252         multTPC++;
1253     }
1254     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
1255         AliAODTrack* trackAOD = event->GetTrack(iTracks);
1256         if (!trackAOD) continue;
1257         if (!(trackAOD->TestFilterBit(16))) continue;
1258         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;
1259         Double_t b[2] = {-99., -99.};
1260         Double_t bCov[3] = {-99., -99., -99.};
1261         if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
1262         if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
1263         multGlob++;
1264     } //track loop
1265     fMultCorBeforeCuts->Fill(multGlob, multTPC);//before all cuts...even before centrality selectrion
1266     //============================================================================================================================
1267     // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
1268     if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
1269     fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
1270     //   cout << "--------------Centrality evaluated-------------------------"<<endl;
1271     if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
1272     {
1273         fCentralityNoPass->Fill(fCentrality);
1274         //    cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
1275         centralitypass = kFALSE;
1276     }else
1277     {
1278         //    cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
1279         centralitypass = kTRUE;
1280     }
1281     if (centralitypass){
1282         fMultCorAfterCentrBeforeCuts->Fill(multGlob, multTPC);
1283         fCentralityBeforePileup->Fill(fCentrality);
1284     }//...after centrality selectrion
1285     //============================================================================================================================
1286     //to remove the bias introduced by multeplicity outliers---------------------
1287     Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
1288     Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
1289     if (TMath::Abs(centv0 - centTrk) > 5.0){
1290         centralitypass = kFALSE;
1291         fCentralityNoPass->Fill(fCentrality);
1292     }
1293     if (centralitypass){
1294         fMultCorAfterVZTRKComp->Fill(multGlob, multTPC);
1295         fCentralityAfterVZTRK->Fill(fCentrality);
1296     }//...after centrality selectrion
1297     //============================================================================================================================
1298     if(fMultCut){
1299         if(fTrigger==1 || fTrigger==4){
1300             if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
1301                 //   cout <<" Trigger ==" <<fTrigger<< endl;
1302                 centralitypass = kFALSE;
1303                 fCentralityNoPass->Fill(fCentrality);
1304             }//2011 Semicentral
1305         }
1306         if(fTrigger==0){
1307             if(! (multTPC > (77.9 + 1.395*multGlob) && multTPC < (187.3 + 1.665*multGlob))){
1308                 //     cout <<" Trigger ==" <<fTrigger<< endl;
1309                 centralitypass = kFALSE;
1310                 fCentralityNoPass->Fill(fCentrality);
1311             }//2011
1312         }//2011 Central
1313     }
1314     if (centralitypass){
1315         fMultCorAfterCorrCut->Fill(multGlob, multTPC);
1316         fCentralityAfterCorrCut->Fill(fCentrality);
1317     }//...after CORR CUT
1318     //=================================All cuts are passed==================++++==================================================
1319     //=================================Now Centrality flattening for central trigger==================++++==================================================
1320     if(fTrigger==0 || fTrigger==4){
1321         if(!IsEventSelectedForCentrFlattening(fCentrality)){
1322             centralitypass = kFALSE;
1323             fCentralityNoPassForFlattening->Fill(fCentrality);
1324         }
1325     }
1326     //==============================fill histo after all cuts==============================++++==================================================
1327     if(centralitypass){
1328         fCentralityPass->Fill(fCentrality);
1329         fMultCorAfterCuts->Fill(multGlob, multTPC);
1330         fMultvsCentr->Fill(fCentrality, multTPC);
1331     }
1332 }
1333 //_____________________________________________________________________________
1334 void AliAnalysisTaskFlowTPCEMCalQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
1335 {
1336     // Set a centrality range ]min, max] and define the method to use for centrality selection
1337     fCentralityMin = CentralityMin;
1338     fCentralityMax = CentralityMax;
1339     fkCentralityMethod = CentralityMethod;
1340 }
1341 //_____________________________________________________________________________
1342 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)
1343 {
1344     //Set ID cuts
1345     fminTPC = minTPC;
1346     fmaxTPC = maxTPC;
1347     fminEovP = minEovP;
1348     fmaxEovP = maxEovP;
1349     fminM20 = minM20;
1350     fmaxM20 = maxM20;
1351     fminM02 = minM02;
1352     fmaxM02 = maxM02;
1353     fDispersion = Dispersion;
1354 }
1355 //_____________________________________________________________________________
1356 //_____________________________________________________________________________
1357 void AliAnalysisTaskFlowTPCEMCalQCSP::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
1358     // set the histo for centrality flattening
1359     // the centrality is flatten in the range minCentr,maxCentr
1360     // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
1361     //                positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
1362     //                negative, h(bin with max in range)*centrRef is used to define the reference (-> defines the maximum loss of events, also in this case the distribution might be not flat)
1363     // switchTRand is used to set the unerflow bin of the histo: if it is < -1 in the analysis the random event selection will be done on using TRandom
1364
1365     if(maxCentr<minCentr){
1366         AliWarning("AliAnalysisCheckCorrdist::Wrong centralities values while setting the histogram for centrality flattening");
1367     }
1368
1369     if(fHistCentrDistr)delete fHistCentrDistr;
1370     fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
1371     fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
1372     Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
1373     Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
1374     fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
1375     Double_t ref=0.,bincont=0.,binrefwidth=1.;
1376     Int_t binref=0;
1377     if(TMath::Abs(centrRef)<0.0001){
1378         binref=fHistCentrDistr->GetMinimumBin();
1379         binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1380         ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1381     }
1382     else if(centrRef>0.){
1383         binref=h->FindBin(centrRef);
1384         if(binref<1||binref>h->GetNbinsX()){
1385             AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
1386         }
1387         binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1388         ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1389     }
1390     else{
1391         if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
1392         binref=fHistCentrDistr->GetMaximumBin();
1393         binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1394         ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
1395     }
1396
1397     for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
1398         if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
1399             bincont=h->GetBinContent(j);
1400             fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
1401             fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
1402         }
1403         else{
1404             h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
1405         }
1406     }
1407
1408     fHistCentrDistr->SetBinContent(0,switchTRand);
1409     return;
1410
1411 }
1412
1413 //-------------------------------------------------
1414 Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::IsEventSelectedForCentrFlattening(Float_t centvalue){
1415     //
1416     //  Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
1417     //  Can be faster if it was required that fHistCentrDistr covers
1418     //  exactly the desired centrality range (e.g. part of the lines below should be done during the
1419     // setting of the histo) and TH1::SetMinimum called
1420     //
1421
1422     if(!fHistCentrDistr) return kTRUE;
1423     // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
1424     //   if(maxbin>fHistCentrDistr->GetNbinsX()){
1425     //     AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
1426     //   }
1427
1428     Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
1429     Double_t bincont=fHistCentrDistr->GetBinContent(bin);
1430     Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
1431
1432     if(fHistCentrDistr->GetBinContent(0)<-0.9999){
1433         if(gRandom->Uniform(1.)<bincont)return kTRUE;
1434         return kFALSE;
1435     }
1436
1437     if(centDigits*100.<bincont)return kTRUE;
1438     return kFALSE;
1439
1440 }
1441 //---------------------------------------------------------------------------
1442