PWGHF/hfe/AliAnalysisTaskFlowITSTPCTOFQCSP.cxx

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