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(kFALSE);
 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             }
 648         }
 649         //==========================================================================================================
 650         //===============================From here eovP cut is used fro QC, SP and EPV0=============================
 651         if(fEovP < fminEovP || fEovP >fmaxEovP) continue;
 652         //==========================================================================================================
 653         //============================Event Plane Method with V0====================================================
 654         Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
 655         Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
 656         Double_t v2Phi[3] = {
 657             v2PhiV0A,
 658             v2PhiV0C,
 659             pt};
 660         fV2Phi->Fill(v2Phi);
 661
 662         Double_t v2PhiVz = TMath::Cos(2*(phi - evPlAngV0));
 663         Double_t v2PhiV0tot[2] = {
 664             v2PhiVz,
 665             pt};
 666         fV2Phivzerotot->Fill(v2PhiV0tot);
 667         //=========================================================================================================
 668         fTPCnsigmaAft->Fill(p,fTPCnSigma);
 669         fInclusiveElecPt->Fill(pt);
 670         fPhi->Fill(phi);
 671         fEta->Fill(eta);
 672         //----------------------Flow of Inclusive Electrons--------------------------------------------------------
 673         AliFlowTrack *sTrack = new AliFlowTrack();
 674         sTrack->Set(track);
 675         sTrack->SetID(track->GetID());
 676         sTrack->SetForRPSelection(kFALSE);
 677         sTrack->SetForPOISelection(kTRUE);
 678         sTrack->SetMass(263732);
 679         for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
 680         {
 681             //   cout << " no of rps " << iRPs << endl;
 682             AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
 683             if (!iRP) continue;
 684             if (!iRP->InRPSelection()) continue;
 685             if( sTrack->GetID() == iRP->GetID())
 686             {
 687                 if(fDebug) printf(" was in RP set");
 688                 //       cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set" <<endl;
 689                 iRP->SetForRPSelection(kFALSE);
 690                 fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
 691             }
 692         } //end of for loop on RPs
 693         fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
 694
 695
 696
 697         if(fDCA){
 698             //----------------------Selection of Photonic Electrons DCA-----------------------------
 699             fNonHFE = new AliSelectNonHFE();
 700             fNonHFE->SetAODanalysis(kTRUE);
 701             fNonHFE->SetInvariantMassCut(fInvmassCut);
 702             if(fOP_angle) fNonHFE->SetOpeningAngleCut(fOpeningAngleCut);
 703             //fNonHFE->SetChi2OverNDFCut(fChi2Cut);
 704             //if(fDCAcutFlag) fNonHFE->SetDCACut(fDCAcut);
 705             fNonHFE->SetAlgorithm("DCA"); //KF
 706             fNonHFE->SetPIDresponse(pidResponse);
 707             fNonHFE->SetTrackCuts(-3,3);
 708
 709             fNonHFE->SetHistAngleBack(fOpeningAngleLS);
 710             fNonHFE->SetHistAngle(fOpeningAngleULS);
 711             //fNonHFE->SetHistDCABack(fDCABack);
 712             //fNonHFE->SetHistDCA(fDCA);
 713             fNonHFE->SetHistMassBack(fInvmassLS1);
 714             fNonHFE->SetHistMass(fInvmassULS1);
 715
 716             fNonHFE->FindNonHFE(iTracks,track,fAOD);
 717
 718             // Int_t *fUlsPartner = fNonHFE->GetPartnersULS();
 719             // Int_t *fLsPartner = fNonHFE->GetPartnersLS();
 720             // Bool_t fUlsIsPartner = kFALSE;
 721             // Bool_t fLsIsPartner = kFALSE;
 722             if(fNonHFE->IsULS()){
 723                 for(Int_t kULS =0; kULS < fNonHFE->GetNULS(); kULS++){
 724                     fULSElecPt->Fill(track->Pt());
 725                 }
 726             }
 727
 728             if(fNonHFE->IsLS()){
 729                 for(Int_t kLS =0; kLS < fNonHFE->GetNLS(); kLS++){
 730                     fLSElecPt->Fill(track->Pt());
 731                 }
 732             }
 733         }
 734
 735         if(!fDCA){
 736             //----------------------Selection of Photonic Electrons KFParticle-----------------------------
 737             Bool_t fFlagPhotonicElec = kFALSE;
 738             SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
 739             if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
 740             // Semi inclusive electron
 741             if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
 742         }
 743
 744
 745
 746     }//end loop on track
 747
 748     PostData(1, fOutputList);
 749     PostData(2, fFlowEvent);
 750     if(fSideBandsFlow){
 751         PostData(3, fFlowEventCont);
 752     }
 753
 754     //----------hfe end---------
 755 }
 756 //_________________________________________
 757 void AliAnalysisTaskFlowTPCEMCalQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track, Bool_t &fFlagPhotonicElec)
 758 {
 759     //Identify non-heavy flavour electrons using Invariant mass method KF
 760
 761     Bool_t flagPhotonicElec = kFALSE;
 762
 763     for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
 764         AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
 765         if (!trackAsso) {
 766             printf("ERROR: Could not receive track %d\n", jTracks);
 767             continue;
 768         }
 769         //  if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
 770         if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
 771         //    if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit) || (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
 772
 773         if(fAssoITSRefit){
 774             if(!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)) continue;
 775         }
 776
 777         if(!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)) continue;
 778
 779         if(jTracks == itrack) continue;
 780         Double_t ptAsso=-999., nsigma=-999.0;
 781         Double_t mass=-999., width = -999;
 782         Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 783         Double_t openingAngle = -999.;
 784         Double_t ptcutonmasshighpt = track->Pt();
 785
 786         ptAsso = trackAsso->Pt();
 787         Short_t chargeAsso = trackAsso->Charge();
 788         Short_t charge = track->Charge();
 789
 790         nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
 791
 792         //80
 793         if(trackAsso->GetTPCNcls() < fAssoTPCCluster) continue;
 794         if(nsigma < -3 || nsigma > 3) continue;
 795         if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
 796         if(ptAsso <0.3) continue;
 797
 798         Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 799         if(charge>0) fPDGe1 = -11;
 800         if(chargeAsso>0) fPDGe2 = -11;
 801
 802         if(charge == chargeAsso) fFlagLS = kTRUE;
 803         if(charge != chargeAsso) fFlagULS = kTRUE;
 804
 805         AliKFParticle::SetField(fAOD->GetMagneticField());
 806         AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
 807         AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
 808         AliKFParticle recg(ge1, ge2);
 809
 810         if(recg.GetNDF()<1) continue;
 811         Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 812         if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 813         recg.GetMass(mass,width);
 814
 815         openingAngle = ge1.GetAngle(ge2);
 816         if(fFlagLS) fOpeningAngleLS->Fill(openingAngle);
 817         if(fFlagULS) fOpeningAngleULS->Fill(openingAngle);
 818         if(fOP_angle)if(openingAngle > fOpeningAngleCut) continue;
 819
 820
 821         if(fFlagLS) fInvmassLS1->Fill(mass);
 822         if(fFlagULS) fInvmassULS1->Fill(mass);
 823
 824         if(ptcutonmasshighpt >= 8.){
 825             if(fFlagLS) fInvmassLS1highpt->Fill(mass);
 826             if(fFlagULS) fInvmassULS1highpt->Fill(mass);
 827         }
 828
 829
 830         if(mass<fInvmassCut){
 831             if(fFlagULS){fULSElecPt->Fill(track->Pt());}
 832             if(fFlagLS){fLSElecPt->Fill(track->Pt());}
 833         }
 834
 835         if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
 836             flagPhotonicElec = kTRUE;
 837         }
 838     }//track loop
 839     fFlagPhotonicElec = flagPhotonicElec;
 840 }
 841 //__________________________________________________________________________________
 842 void AliAnalysisTaskFlowTPCEMCalQCSP::UserCreateOutputObjects()
 843 {
 844     //Create histograms
 845
 846     //----------hfe initialising begin---------
 847     fNullCuts = new AliFlowTrackCuts("null_cuts");
 848
 849     AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 850     cc->SetNbinsMult(10000);
 851     cc->SetMultMin(0);
 852     cc->SetMultMax(10000);
 853
 854     cc->SetNbinsPt(100);
 855     cc->SetPtMin(0);
 856     cc->SetPtMax(50);
 857
 858     cc->SetNbinsPhi(180);
 859     cc->SetPhiMin(0.0);
 860     cc->SetPhiMax(TMath::TwoPi());
 861
 862     cc->SetNbinsEta(30);
 863     cc->SetEtaMin(-7.0);
 864     cc->SetEtaMax(+7.0);
 865
 866     cc->SetNbinsQ(500);
 867     cc->SetQMin(0.0);
 868     cc->SetQMax(3.0);
 869
 870     //--------Initialize PID
 871     fPID->SetHasMCData(kFALSE);
 872     if(!fPID->GetNumberOfPIDdetectors())
 873     {
 874         fPID->AddDetector("TPC", 0);
 875         fPID->AddDetector("EMCAL", 1);
 876     }
 877
 878     fPID->SortDetectors();
 879     fPIDqa = new AliHFEpidQAmanager();
 880     fPIDqa->Initialize(fPID);
 881
 882     //--------Initialize correction Framework and Cuts
 883     fCFM = new AliCFManager;
 884     const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 885     fCFM->SetNStepParticle(kNcutSteps);
 886     for(Int_t istep = 0; istep < kNcutSteps; istep++)
 887         fCFM->SetParticleCutsList(istep, NULL);
 888
 889     if(!fCuts){
 890         AliWarning("Cuts not available. Default cuts will be used");
 891         fCuts = new AliHFEcuts;
 892         fCuts->CreateStandardCuts();
 893     }
 894
 895     fCuts->SetAOD();
 896     fCuts->Initialize(fCFM);
 897     //----------hfe initialising end--------
 898     //---------Output Tlist
 899     fOutputList = new TList();
 900     fOutputList->SetOwner();
 901     fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 902
 903     fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
 904     fOutputList->Add(fNoEvents);
 905
 906     fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",1000,0,50,200,-10,10);
 907     fOutputList->Add(fTPCnsigma);
 908
 909     fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",1000,0,50,200,-10,10);
 910     fOutputList->Add(fTPCnsigmaAft);
 911
 912     fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",1000,0,50,100,0,2);
 913     fOutputList->Add(fTrkEovPBef);
 914
 915     //    fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",1000,0,50,150,0,150);
 916     //     fOutputList->Add(fdEdxBef);
 917
 918     fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",1000,0,100);
 919     fOutputList->Add(fInclusiveElecPt);
 920
 921     fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",1000,0,100);
 922     fOutputList->Add(fPhotoElecPt);
 923
 924     fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",1000,0,100);
 925     fOutputList->Add(fSemiInclElecPt);
 926
 927     fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",1000,0,100);
 928     fOutputList->Add(fULSElecPt);
 929
 930     fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",1000,0,100);
 931     fOutputList->Add(fLSElecPt);
 932
 933     fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
 934     fOutputList->Add(fInvmassLS1);
 935
 936     fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
 937     fOutputList->Add(fInvmassULS1);
 938
 939     fInvmassLS1highpt = new TH1F("fInvmassLS1highpt", "Inv mass of LS (e,e); mass(GeV/c^2) highpt; counts;", 1000,0,1.0);
 940     fOutputList->Add(fInvmassLS1highpt);
 941
 942     fInvmassULS1highpt = new TH1F("fInvmassULS1highpt", "Inv mass of ULS (e,e); mass(GeV/c^2) highpt; counts;", 1000,0,1.0);
 943     fOutputList->Add(fInvmassULS1highpt);
 944
 945     fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
 946     fOutputList->Add(fCentralityPass);
 947
 948     fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
 949     fOutputList->Add(fCentralityNoPass);
 950
 951     fCentralityNoPassForFlattening = new TH1F("fCentralityNoPassForFlattening", "Centrality No Pass for flattening", 101, -1, 100);
 952     fOutputList->Add(fCentralityNoPassForFlattening);
 953
 954     fCentralityBeforePileup = new TH1F("fCentralityBeforePileup", "fCentralityBeforePileup Pass", 101, -1, 100);
 955     fOutputList->Add(fCentralityBeforePileup);
 956
 957     fCentralityAfterVZTRK = new TH1F("fCentralityAfterVZTRK", "fCentralityAfterVZTRK Pass", 101, -1, 100);
 958     fOutputList->Add(fCentralityAfterVZTRK);
 959
 960     fCentralityAfterCorrCut = new TH1F("fCentralityAfterCorrCut", "fCentralityAfterCorrCut Pass", 101, -1, 100);
 961     fOutputList->Add(fCentralityAfterCorrCut);
 962
 963     fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
 964     fOutputList->Add(fPhi);
 965
 966     fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
 967     fOutputList->Add(fEta);
 968
 969     fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
 970     fOutputList->Add(fVZEROA);
 971
 972     fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
 973     fOutputList->Add(fVZEROC);
 974
 975     fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
 976     fOutputList->Add(fTPCM);
 977
 978     fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
 979     fOutputList->Add(fvertex);
 980
 981     fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
 982     fOutputList->Add(fMultCorBeforeCuts);
 983
 984     fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
 985     fOutputList->Add(fMultCorAfterCuts);
 986
 987     fMultCorAfterCentrBeforeCuts = new TH2F("fMultCorAfterCentrBeforeCuts", "TPC vs Global multiplicity (After CC before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
 988     fOutputList->Add(fMultCorAfterCentrBeforeCuts);
 989
 990     fMultCorAfterVZTRKComp = new TH2F("fMultCorAfterVZTRKComp", "TPC vs Global multiplicity (After V0-TRK); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
 991     fOutputList->Add(fMultCorAfterVZTRKComp);
 992
 993     fMultCorAfterCorrCut = new TH2F("fMultCorAfterCorrCut", "TPC vs Global multiplicity (After CorrCut); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
 994     fOutputList->Add(fMultCorAfterCorrCut);
 995
 996     fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
 997     fOutputList->Add(fMultvsCentr);
 998
 999     fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1);
1000     fOutputList->Add(fOpeningAngleLS);
1001
1002     fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1);
1003     fOutputList->Add(fOpeningAngleULS);
1004
1005
1006     //----------------------------------------------------------------------------
1007     EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
1008     fOutputList->Add(EPVzA);
1009     EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
1010     fOutputList->Add(EPVzC);
1011     EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
1012     fOutputList->Add(EPTPC);
1013
1014
1015     EPVz = new TH1D("EPVz", "EPVz", 80, -2, 2);
1016     fOutputList->Add(EPVz);
1017     EPTPCp = new TH1D("EPTPCp", "EPTPCp", 80, -2, 2);
1018     fOutputList->Add(EPTPCp);
1019     EPTPCn = new TH1D("EPTPCn", "EPTPCn", 80, -2, 2);
1020     fOutputList->Add(EPTPCn);
1021
1022     //----------------------------------------------------------------------------
1023     fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
1024     fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1025     fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1026     fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1027     fOutputList->Add(fSubEventDPhiv2);
1028
1029
1030
1031     fSubEventDPhiv2new = new TProfile("fSubEventDPhiv2new", "fSubEventDPhiv2new", 3, 0, 3);
1032     fSubEventDPhiv2new->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
1033     fSubEventDPhiv2new->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
1034     fSubEventDPhiv2new->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
1035     fOutputList->Add(fSubEventDPhiv2new);
1036     //================================Event Plane with VZERO A & C=====================
1037     const Int_t nPtBins = 12;
1038     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};
1039     // v2A, v2C, pt
1040     Int_t    bins[3] = {  50,  50, nPtBins};
1041     Double_t xmin[3] = { -1., -1.,   0};
1042     Double_t xmax[3] = {  1.,  1.,   13};
1043     fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
1044     // Set bin limits for axes which are not standard binned
1045     fV2Phi->SetBinEdges(2, binsPt);
1046     // set axes titles
1047     fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
1048     fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
1049     fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1050     fOutputList->Add(fV2Phi);
1051
1052
1053
1054     //================================Event Plane with VZERO=====================
1055     // const Int_t nPtBins = 10;
1056     // 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};
1057     // v2, pt
1058     Int_t    binsV[2] = {  50,  nPtBins};
1059     Double_t xminV[2] = { -1.,   0};
1060     Double_t xmaxV[2] = {  1.,   13};
1061     fV2Phivzerotot = new THnSparseF("fV2Phivzerotot", "v2:pt", 2, binsV, xminV, xmaxV);
1062     // Set bin limits for axes which are not standard binned
1063     fV2Phivzerotot->SetBinEdges(1, binsPt);
1064     // set axes titles
1065     fV2Phivzerotot->GetAxis(0)->SetTitle("v_{2} (V0)");
1066     fV2Phivzerotot->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
1067     fOutputList->Add(fV2Phivzerotot);
1068
1069
1070
1071     //----------------------------------------------------------------------------
1072
1073     if(fPhiminusPsi){
1074         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)
1075         Double_t xminvElectH[8]={0,    0, -10 ,  0,    0,    0,    0,           0};
1076         Double_t xmaxvElectH[8]={20,   2,  10 ,  2,    2,    2,  TMath::Pi(), TMath::Pi()};
1077         fSparseElectronHadron = new THnSparseD("ElectronHadron","ElectronHadron",8,binsvElectH,xminvElectH,xmaxvElectH);
1078         fSparseElectronHadron->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1079         fSparseElectronHadron->GetAxis(1)->SetTitle("EovP");
1080         fSparseElectronHadron->GetAxis(2)->SetTitle("TPCnSigma");
1081         fSparseElectronHadron->GetAxis(3)->SetTitle("M20");
1082         fSparseElectronHadron->GetAxis(4)->SetTitle("M02");
1083         fSparseElectronHadron->GetAxis(5)->SetTitle("Disp");
1084         fSparseElectronHadron->GetAxis(6)->SetTitle("phiminuspsi V0A");
1085         fSparseElectronHadron->GetAxis(7)->SetTitle("phiminuspsi V0C");
1086         fOutputList->Add(fSparseElectronHadron);
1087     }
1088     //----------------------------------------------------------------------------
1089     //----------------------------------------------------------------------------
1090     if(fpurity){
1091         Int_t binsvpurity[3]={   600,200, 200}; //pt, E/p,TPCnSigma
1092         Double_t xminvpurity[3]={0,    0, -10};
1093         Double_t xmaxvpurity[3]={30,   2,  10};
1094         fSparseElectronpurity = new THnSparseD("Electronpurity","Electronpurity",3,binsvpurity,xminvpurity,xmaxvpurity);
1095         fSparseElectronpurity->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
1096         fSparseElectronpurity->GetAxis(1)->SetTitle("EovP");
1097         fSparseElectronpurity->GetAxis(2)->SetTitle("TPCnSigma");
1098         fOutputList->Add(fSparseElectronpurity);
1099     }
1100     //----------------------------------------------------------------------------
1101
1102     PostData(1,fOutputList);
1103     // create and post flowevent
1104     fFlowEvent = new AliFlowEvent(10000);
1105     PostData(2, fFlowEvent);
1106
1107     if(fSideBandsFlow){
1108         fFlowEventCont = new AliFlowEvent(10000);
1109         PostData(3, fFlowEventCont);
1110     }
1111 }
1112
1113 //________________________________________________________________________
1114 void AliAnalysisTaskFlowTPCEMCalQCSP::Terminate(Option_t *)
1115 {
1116     // Info("Terminate");
1117     AliAnalysisTaskSE::Terminate();
1118 }
1119 //_____________________________________________________________________________
1120 template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::PlotVZeroMultiplcities(const T* event) const
1121 {
1122     // QA multiplicity plots
1123     fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
1124     fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
1125 }
1126 //_____________________________________________________________________________
1127 template <typename T> void AliAnalysisTaskFlowTPCEMCalQCSP::SetNullCuts(T* event)
1128 {
1129     //Set null cuts
1130     if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
1131     fCutsRP->SetEvent(event, MCEvent());
1132     fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
1133     fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
1134     fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
1135     fNullCuts->SetEvent(event, MCEvent());
1136 }
1137 //_____________________________________________________________________________
1138 void AliAnalysisTaskFlowTPCEMCalQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
1139 {
1140     //Prepare flow events
1141     FlowEv->ClearFast();
1142     FlowEv->Fill(fCutsRP, fNullCuts);
1143     FlowEv->SetReferenceMultiplicity(iMulti);
1144     FlowEv->DefineDeadZone(0, 0, 0, 0);
1145     //  FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
1146 }
1147 //_____________________________________________________________________________
1148 Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
1149 {
1150     // Check single track cuts for a given cut step
1151     const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1152     if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1153     return kTRUE;
1154 }
1155 //_________________________________________
1156 void AliAnalysisTaskFlowTPCEMCalQCSP::CheckCentrality(AliAODEvent* event, Bool_t &centralitypass)
1157 {
1158     //============================Multiplicity TPV vs Global===============================================================================
1159     const Int_t nGoodTracks = event->GetNumberOfTracks();
1160     Float_t multTPC(0.); // tpc mult estimate
1161     Float_t multGlob(0.); // global multiplicity
1162     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
1163         AliAODTrack* trackAOD = event->GetTrack(iTracks);
1164         if (!trackAOD) continue;
1165         if (!(trackAOD->TestFilterBit(1))) continue;
1166         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;
1167         multTPC++;
1168     }
1169     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
1170         AliAODTrack* trackAOD = event->GetTrack(iTracks);
1171         if (!trackAOD) continue;
1172         if (!(trackAOD->TestFilterBit(16))) continue;
1173         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;
1174         Double_t b[2] = {-99., -99.};
1175         Double_t bCov[3] = {-99., -99., -99.};
1176         if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
1177         if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
1178         multGlob++;
1179     } //track loop
1180     fMultCorBeforeCuts->Fill(multGlob, multTPC);//before all cuts...even before centrality selectrion
1181     //============================================================================================================================
1182     // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
1183     if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
1184     fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
1185     //   cout << "--------------Centrality evaluated-------------------------"<<endl;
1186     if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
1187     {
1188         fCentralityNoPass->Fill(fCentrality);
1189         //    cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
1190         centralitypass = kFALSE;
1191     }else
1192     {
1193         //    cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
1194         centralitypass = kTRUE;
1195     }
1196     if (centralitypass){
1197         fMultCorAfterCentrBeforeCuts->Fill(multGlob, multTPC);
1198         fCentralityBeforePileup->Fill(fCentrality);
1199     }//...after centrality selectrion
1200     //============================================================================================================================
1201     //to remove the bias introduced by multeplicity outliers---------------------
1202     Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
1203     Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
1204     if (TMath::Abs(centv0 - centTrk) > 5.0){
1205         centralitypass = kFALSE;
1206         fCentralityNoPass->Fill(fCentrality);
1207     }
1208     if (centralitypass){
1209         fMultCorAfterVZTRKComp->Fill(multGlob, multTPC);
1210         fCentralityAfterVZTRK->Fill(fCentrality);
1211     }//...after centrality selectrion
1212     //============================================================================================================================
1213     if(fMultCut){
1214         if(fTrigger==1 || fTrigger==4){
1215             if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
1216                 //   cout <<" Trigger ==" <<fTrigger<< endl;
1217                 centralitypass = kFALSE;
1218                 fCentralityNoPass->Fill(fCentrality);
1219             }//2011 Semicentral
1220         }
1221         if(fTrigger==0){
1222             if(! (multTPC > (77.9 + 1.395*multGlob) && multTPC < (187.3 + 1.665*multGlob))){
1223                 //     cout <<" Trigger ==" <<fTrigger<< endl;
1224                 centralitypass = kFALSE;
1225                 fCentralityNoPass->Fill(fCentrality);
1226             }//2011
1227         }//2011 Central
1228     }
1229     if (centralitypass){
1230         fMultCorAfterCorrCut->Fill(multGlob, multTPC);
1231         fCentralityAfterCorrCut->Fill(fCentrality);
1232     }//...after CORR CUT
1233     //=================================All cuts are passed==================++++==================================================
1234     //=================================Now Centrality flattening for central trigger==================++++==================================================
1235     if(fTrigger==0 || fTrigger==4){
1236         if(!IsEventSelectedForCentrFlattening(fCentrality)){
1237             centralitypass = kFALSE;
1238             fCentralityNoPassForFlattening->Fill(fCentrality);
1239         }
1240     }
1241     //==============================fill histo after all cuts==============================++++==================================================
1242     if(centralitypass){
1243         fCentralityPass->Fill(fCentrality);
1244         fMultCorAfterCuts->Fill(multGlob, multTPC);
1245         fMultvsCentr->Fill(fCentrality, multTPC);
1246     }
1247 }
1248 //_____________________________________________________________________________
1249 void AliAnalysisTaskFlowTPCEMCalQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
1250 {
1251     // Set a centrality range ]min, max] and define the method to use for centrality selection
1252     fCentralityMin = CentralityMin;
1253     fCentralityMax = CentralityMax;
1254     fkCentralityMethod = CentralityMethod;
1255 }
1256 //_____________________________________________________________________________
1257 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)
1258 {
1259     //Set ID cuts
1260     fminTPC = minTPC;
1261     fmaxTPC = maxTPC;
1262     fminEovP = minEovP;
1263     fmaxEovP = maxEovP;
1264     fminM20 = minM20;
1265     fmaxM20 = maxM20;
1266     fminM02 = minM02;
1267     fmaxM02 = maxM02;
1268     fDispersion = Dispersion;
1269 }
1270 //_____________________________________________________________________________
1271 //_____________________________________________________________________________
1272 void AliAnalysisTaskFlowTPCEMCalQCSP::SetHistoForCentralityFlattening(TH1F *h,Double_t minCentr,Double_t maxCentr,Double_t centrRef,Int_t switchTRand){
1273     // set the histo for centrality flattening
1274     // the centrality is flatten in the range minCentr,maxCentr
1275     // if centrRef is zero, the minimum in h within (minCentr,maxCentr) defines the reference
1276     //                positive, the value of h(centrRef) defines the reference (-> the centrality distribution might be not flat in the whole desired range)
1277     //                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)
1278     // 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
1279
1280     if(maxCentr<minCentr){
1281         AliWarning("AliAnalysisCheckCorrdist::Wrong centralities values while setting the histogram for centrality flattening");
1282     }
1283
1284     if(fHistCentrDistr)delete fHistCentrDistr;
1285     fHistCentrDistr=(TH1F*)h->Clone("hCentralityFlat");
1286     fHistCentrDistr->SetTitle("Reference histo for centrality flattening");
1287     Int_t minbin=fHistCentrDistr->FindBin(minCentr*1.00001); // fast if fix bin width
1288     Int_t maxbin=fHistCentrDistr->FindBin(maxCentr*0.9999);
1289     fHistCentrDistr->GetXaxis()->SetRange(minbin,maxbin);
1290     Double_t ref=0.,bincont=0.,binrefwidth=1.;
1291     Int_t binref=0;
1292     if(TMath::Abs(centrRef)<0.0001){
1293         binref=fHistCentrDistr->GetMinimumBin();
1294         binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1295         ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1296     }
1297     else if(centrRef>0.){
1298         binref=h->FindBin(centrRef);
1299         if(binref<1||binref>h->GetNbinsX()){
1300             AliWarning("AliRDHFCuts::Wrong centrality reference value while setting the histogram for centrality flattening");
1301         }
1302         binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1303         ref=fHistCentrDistr->GetBinContent(binref)/binrefwidth;
1304     }
1305     else{
1306         if(centrRef<-1) AliWarning("AliRDHFCuts: with this centrality reference no flattening will be applied");
1307         binref=fHistCentrDistr->GetMaximumBin();
1308         binrefwidth=fHistCentrDistr->GetBinWidth(binref);
1309         ref=fHistCentrDistr->GetMaximum()*TMath::Abs(centrRef)/binrefwidth;
1310     }
1311
1312     for(Int_t j=1;j<=h->GetNbinsX();j++){// Now set the "probabilities"
1313         if(h->GetBinLowEdge(j)*1.0001>=minCentr&&h->GetBinLowEdge(j+1)*0.9999<=maxCentr){
1314             bincont=h->GetBinContent(j);
1315             fHistCentrDistr->SetBinContent(j,ref/bincont*h->GetBinWidth(j));
1316             fHistCentrDistr->SetBinError(j,h->GetBinError(j)*ref/bincont);
1317         }
1318         else{
1319             h->SetBinContent(j,1.1);// prob > 1 to assure that events will not be rejected
1320         }
1321     }
1322
1323     fHistCentrDistr->SetBinContent(0,switchTRand);
1324     return;
1325
1326 }
1327
1328 //-------------------------------------------------
1329 Bool_t AliAnalysisTaskFlowTPCEMCalQCSP::IsEventSelectedForCentrFlattening(Float_t centvalue){
1330     //
1331     //  Random event selection, based on fHistCentrDistr, to flatten the centrality distribution
1332     //  Can be faster if it was required that fHistCentrDistr covers
1333     //  exactly the desired centrality range (e.g. part of the lines below should be done during the
1334     // setting of the histo) and TH1::SetMinimum called
1335     //
1336
1337     if(!fHistCentrDistr) return kTRUE;
1338     // Int_t maxbin=fHistCentrDistr->FindBin(fMaxCentrality*0.9999);
1339     //   if(maxbin>fHistCentrDistr->GetNbinsX()){
1340     //     AliWarning("AliRDHFCuts: The maximum centrality exceeds the x-axis limit of the histogram for centrality flattening");
1341     //   }
1342
1343     Int_t bin=fHistCentrDistr->FindBin(centvalue); // Fast if the histo has a fix bin
1344     Double_t bincont=fHistCentrDistr->GetBinContent(bin);
1345     Double_t centDigits=centvalue-(Int_t)(centvalue*100.)/100.;// this is to extract a random number between 0 and 0.01
1346
1347     if(fHistCentrDistr->GetBinContent(0)<-0.9999){
1348         if(gRandom->Uniform(1.)<bincont)return kTRUE;
1349         return kFALSE;
1350     }
1351
1352     if(centDigits*100.<bincont)return kTRUE;
1353     return kFALSE;
1354
1355 }
1356 //---------------------------------------------------------------------------
1357