PWGHF/hfe/AliAnalysisTaskFlowITSTPCTOFQCSP.cxx

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