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