PWGHF/hfe/AliAnalysisTaskFlowTPCEMCalQCSP.cxx

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