PWGHF/hfe/AliAnalysisTaskFlowTPCEMCalQCSP.cxx

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