PWGHF/hfe/AliAnalysisTaskFlowTPCEMCalQCSP.cxx

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