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