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