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