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