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