PWGHF/hfe/AliAnalysisTaskFlowTPCTOFQCSP.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 "AliAnalysisTaskFlowTPCTOFQCSP.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 "AliAODPid.h"
  97 #include "AliFlowTrack.h"
  98 #include "AliAnalysisTaskVnV0.h"
  99
 100
 101 class AliFlowTrackCuts;
 102
 103 using namespace std;
 104
 105 ClassImp(AliAnalysisTaskFlowTPCTOFQCSP)
 106 //________________________________________________________________________
 107   AliAnalysisTaskFlowTPCTOFQCSP::AliAnalysisTaskFlowTPCTOFQCSP(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 ,fTrigger(0)
 127 ,fPhi(0)
 128 ,fEta(0)
 129 ,fVZEROA(0)
 130 ,fVZEROC(0)
 131 ,fTPCM(0)
 132 ,fNoEvents(0)
 133 ,fInclusiveElecPt(0)
 134 //,fTPCnsigma(0)
 135 ,fTPCnsigmaAft(0)
 136 //,fTOFns(0)
 137 ,fTOFnsAft(0)
 138 ,fTOFBetaAft(0)
 139 ,fCentralityPass(0)
 140 ,fCentralityNoPass(0)
 141 ,fInvmassLS1(0)
 142 ,fInvmassULS1(0)
 143 ,fPhotoElecPt(0)
 144 ,fSemiInclElecPt(0)
 145 ,fULSElecPt(0)
 146 ,fLSElecPt(0)
 147 ,fMultCorAfterCuts(0)
 148 ,fMultvsCentr(0)
 149 ,fSubEventDPhiv2(0)
 150 ,EPVzA(0)
 151 ,EPVzC(0)
 152 ,EPTPC(0)
 153 ,fV2Phi(0)
 154 ,fvertex(0)
 155 ,fMultCorBeforeCuts(0)
 156 ,fPIDResponse(0)
 157 ,fQAPid(0)
 158 ,fminTPCnsigma(0)
 159 ,fmaxTPCnsigma(3)
 160 ,fQAPIDSparse(kFALSE)
 161 ,fminTOFnSigma(-2)
 162 ,fmaxTOFnSigma(2)
 163 ,fQAPidSparse(0)
 164 ,fTPCS(0)
 165 ,fVz(0)
 166 ,fPhiCut(0)
 167 {
 168   //Named constructor
 169
 170 //  fPID = new AliHFEpid("hfePid");
 171   // Define input and output slots here
 172   // Input slot #0 works with a TChain
 173   DefineInput(0, TChain::Class());
 174   // Output slot #0 id reserved by the base class for AOD
 175   // Output slot #1 writes into a TH1 container
 176   // DefineOutput(1, TH1I::Class());
 177   DefineOutput(1, TList::Class());
 178   DefineOutput(2, AliFlowEventSimple::Class());
 179 }
 180
 181 //________________________________________________________________________
 182 AliAnalysisTaskFlowTPCTOFQCSP::AliAnalysisTaskFlowTPCTOFQCSP()
 183   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisElectFlow")
 184 ,fDebug(0)
 185 ,fAOD(0)
 186 ,fOutputList(0)
 187 ,fCuts(0)
 188 ,fIdentifiedAsOutInz(kFALSE)
 189 ,fPassTheEventCut(kFALSE)
 190 ,fCFM(0)
 191 //,fPID(0)
 192 //,fPIDqa(0)
 193 ,fCutsRP(0)     // track cuts for reference particles
 194 ,fNullCuts(0) // dummy cuts for flow event tracks
 195 ,fFlowEvent(0) //! flow events (one for each inv mass band)
 196 ,fkCentralityMethod(0)
 197 ,fCentrality(0)
 198 ,fCentralityMin(0)
 199 ,fCentralityMax(0)
 200 ,fInvmassCut(0)
 201 ,fTrigger(0)
 202 ,fPhi(0)
 203 ,fEta(0)
 204 ,fVZEROA(0)
 205 ,fVZEROC(0)
 206 ,fTPCM(0)
 207 ,fNoEvents(0)
 208 ,fInclusiveElecPt(0)
 209 //,fTPCnsigma(0)
 210 ,fTPCnsigmaAft(0)
 211 //,fTOFns(0)
 212 ,fTOFnsAft(0)
 213 ,fTOFBetaAft(0)
 214 ,fCentralityPass(0)
 215 ,fCentralityNoPass(0)
 216 ,fInvmassLS1(0)
 217 ,fInvmassULS1(0)
 218 ,fPhotoElecPt(0)
 219 ,fSemiInclElecPt(0)
 220 ,fULSElecPt(0)
 221 ,fLSElecPt(0)
 222 ,fMultCorAfterCuts(0)
 223 ,fMultvsCentr(0)
 224 ,fSubEventDPhiv2(0)
 225 ,EPVzA(0)
 226 ,EPVzC(0)
 227 ,EPTPC(0)
 228 ,fV2Phi(0)
 229 ,fvertex(0)
 230 ,fMultCorBeforeCuts(0)
 231 ,fPIDResponse(0)
 232 ,fQAPid(0)
 233 ,fminTPCnsigma(0)
 234 ,fmaxTPCnsigma(3)
 235 ,fQAPIDSparse(kFALSE)
 236 ,fminTOFnSigma(-2)
 237 ,fmaxTOFnSigma(2)
 238 ,fQAPidSparse(0)
 239 ,fTPCS(0)
 240 ,fVz(0)
 241 ,fPhiCut(0)
 242 {
 243   //Default constructor
 244 //  fPID = new AliHFEpid("hfePid");
 245   // Constructor
 246   // Define input and output slots here
 247   // Input slot #0 works with a TChain
 248   DefineInput(0, TChain::Class());
 249   // Output slot #0 id reserved by the base class for AOD
 250   // Output slot #1 writes into a TH1 container
 251   // DefineOutput(1, TH1I::Class());
 252   DefineOutput(1, TList::Class());
 253   DefineOutput(2, AliFlowEventSimple::Class());
 254     //DefineOutput(3, TTree::Class());
 255 }
 256 //_________________________________________
 257
 258 AliAnalysisTaskFlowTPCTOFQCSP::~AliAnalysisTaskFlowTPCTOFQCSP()
 259 {
 260   //Destructor
 261
 262   delete fOutputList;
 263 //  delete fPID;
 264   delete fPIDResponse;
 265   delete fCFM;
 266 //  delete fPIDqa;
 267   if (fOutputList) delete fOutputList;
 268   if (fFlowEvent) delete fFlowEvent;
 269
 270 }
 271 //_________________________________________
 272
 273 void AliAnalysisTaskFlowTPCTOFQCSP::UserExec(Option_t*)
 274 {
 275   //Main loop
 276   //Called for each event
 277
 278   // create pointer to event
 279
 280   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 281   if (!fAOD)
 282   {
 283     printf("ERROR: fAOD not available\n");
 284     return;
 285   }
 286
 287   if(!fCuts)
 288   {
 289     AliError("HFE cuts not available");
 290     return;
 291   }
 292
 293 //  if(!fPID->IsInitialized())
 294 //  {
 295     // Initialize PID with the given run number
 296 //    AliWarning("PID not initialised, get from Run no");
 297 //    fPID->InitializePID(fAOD->GetRunNumber());
 298 //  }
 299
 300   // cout << "kTrigger   ==   " << fTrigger <<endl;
 301
 302     if(fTrigger==0){
 303 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral) ) return;
 304     }
 305     if(fTrigger==1){
 306 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kSemiCentral))) return;
 307     }
 308     if(fTrigger==2){
 309 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMCEGA) ) return;
 310     }
 311     if(fTrigger==3){
 312 if(!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) ) return;
 313     }
 314
 315
 316 //---------------CENTRALITY AND EVENT SELECTION-----------------------
 317
 318
 319
 320    Int_t fNOtrks =  fAOD->GetNumberOfTracks();
 321     Float_t vtxz = -999;
 322     const AliAODVertex* trkVtx = fAOD->GetPrimaryVertex();
 323     if (!trkVtx || trkVtx->GetNContributors()<=0)return;
 324     TString vtxTtl = trkVtx->GetTitle();
 325     if (!vtxTtl.Contains("VertexerTracks"))return;
 326     const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
 327     if (!spdVtx || spdVtx->GetNContributors()<=0)return;
 328     if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5)return;
 329     vtxz = trkVtx->GetZ();
 330     if(TMath::Abs(vtxz)>fVz)return;
 331
 332 // Event cut
 333     if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) return;
 334     if(fNOtrks<2) return;
 335
 336
 337     Bool_t pass = kFALSE; //to select centrality
 338     CheckCentrality(fAOD,pass);
 339     if(!pass)return;
 340
 341     fvertex->Fill(vtxz);
 342
 343
 344     fNoEvents->Fill(0);
 345     PlotVZeroMultiplcities(fAOD);
 346
 347     SetNullCuts(fAOD);
 348     PrepareFlowEvent(fAOD->GetNumberOfTracks(),fFlowEvent);    //Calculate event plane Qvector and EP resolution for inclusive
 349
 350
 351     fCFM->SetRecEventInfo(fAOD);
 352
 353   // Look for kink mother
 354   Int_t numberofvertices = fAOD->GetNumberOfVertices();
 355   Double_t listofmotherkink[numberofvertices];
 356   Int_t numberofmotherkink = 0;
 357   for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
 358     AliAODVertex *aodvertex = fAOD->GetVertex(ivertex);
 359     if(!aodvertex) continue;
 360     if(aodvertex->GetType()==AliAODVertex::kKink) {
 361       AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
 362       if(!mother) continue;
 363       Int_t idmother = mother->GetID();
 364       listofmotherkink[numberofmotherkink] = idmother;
 365       //printf("ID %d\n",idmother);
 366       numberofmotherkink++;
 367     }
 368   }
 369
 370 //=============================================V0EP from Alex======================================================================
 371     Double_t qxEPa = 0, qyEPa = 0;
 372     Double_t qxEPc = 0, qyEPc = 0;
 373
 374     Double_t evPlAngV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 8, 2, qxEPa, qyEPa);
 375     Double_t evPlAngV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD, 9, 2, qxEPc, qyEPc);
 376
 377
 378     Double_t Qx2 = 0, Qy2 = 0;
 379
 380     for (Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++){
 381
 382         AliAODTrack* aodTrack = fAOD->GetTrack(iT);
 383
 384         if (!aodTrack)
 385             continue;
 386
 387         if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < 70) || (aodTrack->Pt() >= 20.0))
 388             continue;
 389
 390         if (!aodTrack->TestFilterBit(128))
 391             continue;
 392
 393         Qx2 += TMath::Cos(2*aodTrack->Phi());
 394         Qy2 += TMath::Sin(2*aodTrack->Phi());
 395     }
 396
 397     Double_t evPlAngTPC = TMath::ATan2(Qy2, Qx2)/2.;
 398
 399     EPVzA->Fill(evPlAngV0A);
 400     EPVzC->Fill(evPlAngV0C);
 401     EPTPC->Fill(evPlAngTPC);
 402
 403     fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(evPlAngV0A-evPlAngTPC))); // vzeroa - tpc
 404     fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(evPlAngV0A-evPlAngV0C))); // vzeroa - vzeroc
 405     fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(evPlAngV0C-evPlAngTPC))); // tpc - vzeroc
 406 //====================================================================================================================
 407
 408  AliAODTrack *track = NULL;
 409
 410 // Track loop
 411  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)
 412  {
 413    track = fAOD->GetTrack(iTracks);
 414    if (!track)
 415    {
 416      printf("ERROR: Could not receive track %d\n", iTracks);
 417      continue;
 418    }
 419
 420    if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
 421
 422 //--------------------------------------hfe begin-----------------------------------------------------------
 423 //==========================================================================================================
 424 //======================================track cuts==========================================================
 425    if(track->Eta()<-0.8 || track->Eta()>0.8)    continue;    //eta cuts on candidates
 426
 427    // RecKine: ITSTPC cuts
 428    if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 429
 430    // Reject kink mother
 431    Bool_t kinkmotherpass = kTRUE;
 432    for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
 433      if(track->GetID() == listofmotherkink[kinkmother]) {
 434        kinkmotherpass = kFALSE;
 435        continue;
 436      }
 437    }
 438    if(!kinkmotherpass) continue;
 439
 440    // RecPrim
 441  //  if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;  //deleted for DCA absence
 442    // HFEcuts: ITS layers cuts
 443    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 444    // HFE cuts: TPC PID cleanup
 445    if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 446 //==========================================================================================================
 447      Double_t eta = track->Eta();
 448      Double_t phi = track->Phi();
 449      Double_t pt = track->Pt();         //pt track after cuts
 450      if(pt<1.) continue;
 451 //==========================================================================================================
 452 //=========================================PID==============================================================
 453      if(track->GetTPCsignalN() < fTPCS) continue;
 454      Float_t fTPCnSigma = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
 455      Float_t fTOFnSigma = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kElectron);
 456
 457   //   Float_t eDEDX = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE);
 458
 459 //     fTPCnsigma->Fill(track->P(),fTPCnSigma);
 460 //     fTOFns->Fill(track->P(),fTOFnSigma);
 461
 462      Double_t valPidSparse[4] = {
 463          pt,
 464          fTPCnSigma,
 465          fTOFnSigma,
 466          phi
 467         };
 468      fQAPidSparse->Fill(valPidSparse);
 469
 470      /* || track->GetTPCsignal() < 72 || track->GetTPCsignal() > 95*/
 471
 472      if(fTOFnSigma < fminTOFnSigma || fTOFnSigma > fmaxTOFnSigma) continue;
 473      if(fTPCnSigma < fminTPCnsigma || fTPCnSigma > fmaxTPCnsigma) continue; //cuts on nsigma tpc
 474 //==========================================================================================================
 475 //=========================================QA PID SPARSE====================================================
 476      Float_t timeTOF = track->GetTOFsignal();
 477      Double_t intTime[5] = {-99., -99., -99., -99., -99.};
 478      track->GetIntegratedTimes(intTime);
 479      Float_t timeElec = intTime[0];
 480      Float_t intLength = 2.99792458e-2* timeElec;
 481      Double_t beta = 0.1;
 482      if ((intLength > 0) && (timeTOF > 0))
 483          beta = intLength/2.99792458e-2/timeTOF;
 484
 485     if(fQAPIDSparse){
 486      Double_t valPid[4] = {
 487          track->P(),
 488          track->GetTPCsignal(),
 489          beta,
 490          track->Charge()
 491          };
 492      fQAPid->Fill(valPid);
 493      }
 494
 495
 496       if(fPhiCut){
 497       if(phi<1.4 || phi >3.14)continue; //to have same EMCal phi acceptance
 498       }
 499
 500 //==========================================================================================================
 501 //============================Event Plane Method with V0====================================================
 502      Double_t v2PhiV0A = TMath::Cos(2*(phi - evPlAngV0A));
 503      Double_t v2PhiV0C = TMath::Cos(2*(phi - evPlAngV0C));
 504      Double_t v2Phi[3] = {
 505          v2PhiV0A,
 506          v2PhiV0C,
 507          pt};
 508      fV2Phi->Fill(v2Phi);
 509 //=========================================================================================================
 510    fTPCnsigmaAft->Fill(track->P(),fTPCnSigma);
 511    fTOFnsAft->Fill(track->P(),fTOFnSigma);
 512    fTOFBetaAft->Fill(track->P(),beta);
 513    fInclusiveElecPt->Fill(pt);
 514    fPhi->Fill(phi);
 515    fEta->Fill(eta);
 516 //=========================================================================================================
 517 //----------------------Flow of Inclusive Electrons--------------------------------------------------------
 518    AliFlowTrack *sTrack = new AliFlowTrack();
 519      sTrack->Set(track);
 520      sTrack->SetID(track->GetID());
 521      sTrack->SetForRPSelection(kFALSE);
 522      sTrack->SetForPOISelection(kTRUE);
 523      sTrack->SetMass(263732);
 524      for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
 525      {
 526       //   cout << " no of rps " << iRPs << endl;
 527          AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
 528          if (!iRP) continue;
 529          if (!iRP->InRPSelection()) continue;
 530         if( sTrack->GetID() == iRP->GetID())
 531         {
 532           if(fDebug) printf(" was in RP set");
 533           //  cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set====REMOVED" <<endl;
 534             iRP->SetForRPSelection(kFALSE);
 535             fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
 536         }
 537       } //end of for loop on RPs
 538    fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
 539 //=========================================================================================================
 540 //----------------------Selection and Flow of Photonic Electrons-----------------------------
 541    Bool_t fFlagPhotonicElec = kFALSE;
 542    SelectPhotonicElectron(iTracks,track,fFlagPhotonicElec);
 543    if(fFlagPhotonicElec){fPhotoElecPt->Fill(pt);}
 544    // Semi inclusive electron
 545    if(!fFlagPhotonicElec){fSemiInclElecPt->Fill(pt);}
 546 //-------------------------------------------------------------------------------------------
 547
 548  }//end loop on track
 549  PostData(1, fOutputList);
 550  PostData(2, fFlowEvent);
 551  //----------hfe end---------
 552 }
 553 //_________________________________________
 554 void AliAnalysisTaskFlowTPCTOFQCSP::SelectPhotonicElectron(Int_t itrack,const AliAODTrack *track, Bool_t &fFlagPhotonicElec)
 555 {
 556
 557   //Identify non-heavy flavour electrons using Invariant mass method
 558   Bool_t flagPhotonicElec = kFALSE;
 559
 560   for(Int_t jTracks = 0; jTracks<fAOD->GetNumberOfTracks(); jTracks++){
 561     AliAODTrack *trackAsso = fAOD->GetTrack(jTracks);
 562     if (!trackAsso) {
 563       printf("ERROR: Could not receive track %d\n", jTracks);
 564       continue;
 565     }
 566     //  if(!track->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) continue;  // TESTBIT FOR AOD double Counting
 567       if(!trackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly)) continue;
 568       if((!(trackAsso->GetStatus()&AliESDtrack::kITSrefit)|| (!(trackAsso->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
 569
 570
 571     if(jTracks == itrack) continue;
 572     Double_t ptAsso=-999., nsigma=-999.0;
 573     Double_t mass=-999., width = -999;
 574     Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 575
 576
 577     ptAsso = trackAsso->Pt();
 578     Short_t chargeAsso = trackAsso->Charge();
 579     Short_t charge = track->Charge();
 580    // nsigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000;
 581       nsigma = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
 582
 583     if(trackAsso->GetTPCNcls() < 80) continue;
 584     if(nsigma < -3 || nsigma > 3) continue;
 585     if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9) continue;
 586     if(ptAsso <0.3) continue;
 587
 588     Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 589     if(charge>0) fPDGe1 = -11;
 590     if(chargeAsso>0) fPDGe2 = -11;
 591
 592     if(charge == chargeAsso) fFlagLS = kTRUE;
 593     if(charge != chargeAsso) fFlagULS = kTRUE;
 594
 595     AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
 596     AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
 597     AliKFParticle recg(ge1, ge2);
 598
 599     if(recg.GetNDF()<1) continue;
 600     Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 601     if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 602     recg.GetMass(mass,width);
 603
 604     if(fFlagLS) fInvmassLS1->Fill(mass);
 605     if(fFlagULS) fInvmassULS1->Fill(mass);
 606
 607        if(mass<fInvmassCut){
 608        if(fFlagULS){fULSElecPt->Fill(track->Pt());}
 609        if(fFlagLS){fLSElecPt->Fill(track->Pt());}
 610        }
 611
 612     if(mass<fInvmassCut && fFlagULS && !flagPhotonicElec){
 613       flagPhotonicElec = kTRUE;
 614     }
 615   }//track loop
 616   fFlagPhotonicElec = flagPhotonicElec;
 617 }
 618 //___________________________________________
 619 void AliAnalysisTaskFlowTPCTOFQCSP::UserCreateOutputObjects()
 620 {
 621
 622   //Create histograms
 623   //----------hfe initialising begin---------
 624    fNullCuts = new AliFlowTrackCuts("null_cuts");
 625
 626    AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 627   cc->SetNbinsMult(10000);
 628   cc->SetMultMin(0);
 629   cc->SetMultMax(10000);
 630
 631   cc->SetNbinsPt(100);
 632   cc->SetPtMin(0);
 633   cc->SetPtMax(50);
 634
 635   cc->SetNbinsPhi(180);
 636   cc->SetPhiMin(0.0);
 637   cc->SetPhiMax(TMath::TwoPi());
 638
 639   cc->SetNbinsEta(200);
 640   cc->SetEtaMin(-8.0);
 641   cc->SetEtaMax(+8.0);
 642
 643   cc->SetNbinsQ(500);
 644   cc->SetQMin(0.0);
 645   cc->SetQMax(3.0);
 646
 647
 648     AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 649     AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
 650     if (!inputHandler)
 651         AliFatal("Input handler needed");
 652
 653   //pid response object
 654   fPIDResponse=inputHandler->GetPIDResponse();
 655   if (!fPIDResponse)AliError("PIDResponse object was not created");
 656
 657
 658     //--------Initialize PID
 659     //    fPID->SetHasMCData(kFALSE);
 660     //    if(!fPID->GetNumberOfPIDdetectors())
 661     //   {
 662     //       fPID->AddDetector("TPC", 0);
 663     //   fPID->AddDetector("EMCAL", 1);
 664     //   }
 665
 666     //   fPID->SortDetectors();
 667     //   fPIDqa = new AliHFEpidQAmanager();
 668     //   fPIDqa->Initialize(fPID);
 669
 670
 671
 672   //--------Initialize correction Framework and Cuts
 673   fCFM = new AliCFManager;
 674   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 675   fCFM->SetNStepParticle(kNcutSteps);
 676   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 677     fCFM->SetParticleCutsList(istep, NULL);
 678
 679   if(!fCuts){
 680     AliWarning("Cuts not available. Default cuts will be used");
 681     fCuts = new AliHFEcuts;
 682     fCuts->CreateStandardCuts();
 683   }
 684
 685   fCuts->SetAOD();
 686   fCuts->Initialize(fCFM);
 687   //----------hfe initialising end--------
 688   //---------Output Tlist
 689   fOutputList = new TList();
 690   fOutputList->SetOwner();
 691 //  fOutputList->Add(fPIDqa->MakeList("PIDQA"));
 692
 693   fNoEvents = new TH1F("fNoEvents","",1,0,1) ;
 694   fOutputList->Add(fNoEvents);
 695
 696 //  fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma before HFE pid",1000,0,50,200,-10,10);
 697  // fOutputList->Add(fTPCnsigma);
 698
 699   fTPCnsigmaAft = new TH2F("fTPCnsigmaAft", "TPC - n sigma after HFE pid",1000,0,50,200,-10,10);
 700   fOutputList->Add(fTPCnsigmaAft);
 701
 702 //  fTOFns = new TH2F("fTOFns","track TOFnSigma",1000,0,50,200,-10,10);
 703 //  fOutputList->Add(fTOFns);
 704
 705   fTOFnsAft = new TH2F("fTOFnsAft","track TOFnSigma",1000,0,50,200,-10,10);
 706   fOutputList->Add(fTOFnsAft);
 707
 708   fTOFBetaAft = new TH2F("fTOFBetaAft","track TOFBeta",1000,0,50,120,0,1.2);
 709   fOutputList->Add(fTOFBetaAft);
 710
 711   fInclusiveElecPt = new TH1F("fInclElecPt", "Inclusive electron pt",1000,0,100);
 712   fOutputList->Add(fInclusiveElecPt);
 713
 714   fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",1000,0,100);
 715   fOutputList->Add(fPhotoElecPt);
 716
 717   fSemiInclElecPt = new TH1F("fSemiInclElecPt", "Semi-inclusive electron pt",1000,0,100);
 718   fOutputList->Add(fSemiInclElecPt);
 719
 720   fULSElecPt = new TH1F("fULSElecPt", "ULS electron pt",1000,0,100);
 721   fOutputList->Add(fULSElecPt);
 722
 723   fLSElecPt = new TH1F("fLSElecPt", "LS electron pt",1000,0,100);
 724   fOutputList->Add(fLSElecPt);
 725
 726   fInvmassLS1 = new TH1F("fInvmassLS1", "Inv mass of LS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
 727   fOutputList->Add(fInvmassLS1);
 728
 729   fInvmassULS1 = new TH1F("fInvmassULS1", "Inv mass of ULS (e,e); mass(GeV/c^2); counts;", 1000,0,1.0);
 730   fOutputList->Add(fInvmassULS1);
 731
 732   fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
 733   fOutputList->Add(fCentralityPass);
 734
 735   fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
 736   fOutputList->Add(fCentralityNoPass);
 737
 738   fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
 739   fOutputList->Add(fPhi);
 740
 741   fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
 742   fOutputList->Add(fEta);
 743
 744   fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
 745   fOutputList->Add(fVZEROA);
 746
 747   fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
 748   fOutputList->Add(fVZEROC);
 749
 750   fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
 751   fOutputList->Add(fTPCM);
 752
 753   fvertex = new TH1D("fvertex", "vertex distribution", 300, -15,15);
 754   fOutputList->Add(fvertex);
 755
 756   fMultCorBeforeCuts = new TH2F("fMultCorBeforeCuts", "TPC vs Global multiplicity (Before cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
 757   fOutputList->Add(fMultCorBeforeCuts);
 758
 759   fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
 760   fOutputList->Add(fMultCorAfterCuts);
 761
 762   fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 100, 0., 100, 100, 0, 3000);
 763   fOutputList->Add(fMultvsCentr);
 764
 765 //----------------------------------------------------------------------------
 766     EPVzA = new TH1D("EPVzA", "EPVzA", 80, -2, 2);
 767     fOutputList->Add(EPVzA);
 768     EPVzC = new TH1D("EPVzC", "EPVzC", 80, -2, 2);
 769     fOutputList->Add(EPVzC);
 770     EPTPC = new TH1D("EPTPC", "EPTPC", 80, -2, 2);
 771     fOutputList->Add(EPTPC);
 772 //----------------------------------------------------------------------------
 773     fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
 774     fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
 775     fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
 776     fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
 777     fOutputList->Add(fSubEventDPhiv2);
 778 //================================Event Plane with VZERO=====================
 779     const Int_t nPtBins = 10;
 780     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};
 781     // v2A, v2C, pt
 782     Int_t    bins[3] = {  50,  50, nPtBins};
 783     Double_t xmin[3] = { -1., -1.,   0};
 784     Double_t xmax[3] = {  1.,  1.,   8};
 785     fV2Phi = new THnSparseF("fV2Phi", "v2A:v2C:pt", 3, bins, xmin, xmax);
 786     // Set bin limits for axes which are not standard binned
 787     fV2Phi->SetBinEdges(2, binsPt);
 788     // set axes titles
 789     fV2Phi->GetAxis(0)->SetTitle("v_{2} (V0A)");
 790     fV2Phi->GetAxis(1)->SetTitle("v_{2} (V0C)");
 791     fV2Phi->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 792     fOutputList->Add(fV2Phi);
 793 //----------------------------------------------------------------------------
 794 //----------------------------------------------------------------------------
 795     if(fQAPIDSparse){
 796     Int_t    binsQA[4] = { 150,  100,  120,    3};
 797     Double_t xminQA[4] = { 0.,    50,   0, -1.5};
 798     Double_t xmaxQA[4] = { 15.,  150, 1.2,  1.5};
 799     fQAPid = new THnSparseF("fQAPid", "p:dEdx:beta:ch", 4, binsQA, xminQA, xmaxQA);
 800     fQAPid->GetAxis(0)->SetTitle("p (Gev/c");
 801     fQAPid->GetAxis(1)->SetTitle("dE/dx");
 802     fQAPid->GetAxis(2)->SetTitle("#beta (TOF)");
 803     fQAPid->GetAxis(3)->SetTitle("charge");
 804     fOutputList->Add(fQAPid);
 805     }
 806 //===========================================================================
 807         Int_t    binsQA2[5] = { 200,  200, 200, 100};
 808         Double_t xminQA2[5] = { 0.,   -10, -4, 0};
 809         Double_t xmaxQA2[5] = { 20.,   10,  4, 6.28};
 810         fQAPidSparse = new THnSparseF("fQAPidSparse", "pt:p:tpcnsigma:tofnsigma:phi", 4, binsQA2, xminQA2, xmaxQA2);
 811         fQAPidSparse->GetAxis(0)->SetTitle("pt (Gev/c");
 812         fQAPidSparse->GetAxis(1)->SetTitle("tpc nsigma");
 813         fQAPidSparse->GetAxis(2)->SetTitle("tofnsigma");
 814         fQAPidSparse->GetAxis(3)->SetTitle("phi");
 815         fOutputList->Add(fQAPidSparse);
 816 //===========================================================================
 817   PostData(1,fOutputList);
 818  // create and post flowevent
 819   fFlowEvent = new AliFlowEvent(10000);
 820   PostData(2, fFlowEvent);
 821  }
 822 //________________________________________________________________________
 823 void AliAnalysisTaskFlowTPCTOFQCSP::Terminate(Option_t *)
 824 {
 825   // Info("Terminate");
 826   AliAnalysisTaskSE::Terminate();
 827 }
 828 //_____________________________________________________________________________
 829 template <typename T> void AliAnalysisTaskFlowTPCTOFQCSP::PlotVZeroMultiplcities(const T* event) const
 830 {
 831   // QA multiplicity plots
 832   fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
 833   fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
 834 }
 835 //_____________________________________________________________________________
 836 template <typename T> void AliAnalysisTaskFlowTPCTOFQCSP::SetNullCuts(T* event)
 837 {
 838  //Set null cuts
 839     if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
 840     fCutsRP->SetEvent(event, MCEvent());
 841     fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
 842     fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
 843     fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
 844     fNullCuts->SetEvent(event, MCEvent());
 845 }
 846 //_____________________________________________________________________________
 847 void AliAnalysisTaskFlowTPCTOFQCSP::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
 848 {
 849  //Prepare flow events
 850    FlowEv->ClearFast();
 851    FlowEv->Fill(fCutsRP, fNullCuts);
 852    FlowEv->SetReferenceMultiplicity(iMulti);
 853    FlowEv->DefineDeadZone(0, 0, 0, 0);
 854  //  FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
 855 }
 856 //_____________________________________________________________________________
 857 Bool_t AliAnalysisTaskFlowTPCTOFQCSP::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 858 {
 859   // Check single track cuts for a given cut step
 860   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 861   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 862   return kTRUE;
 863 }
 864 //_________________________________________
 865 void AliAnalysisTaskFlowTPCTOFQCSP::CheckCentrality(AliAODEvent* event, Bool_t &centralitypass)
 866 {
 867   // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
 868   if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
 869   fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
 870  // cout << "--------------Centrality evaluated-------------------------"<<endl;
 871   if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax))
 872   {
 873     fCentralityNoPass->Fill(fCentrality);
 874     //cout << "--------------Fill no pass-----"<< fCentrality <<"--------------------"<<endl;
 875     centralitypass = kFALSE;
 876   }else
 877   {
 878   //  cout << "--------------Fill pass----"<< fCentrality <<"---------------------"<<endl;
 879     centralitypass = kTRUE;
 880   }
 881 //to remove the bias introduced by multeplicity outliers---------------------
 882     Float_t centTrk = event->GetCentrality()->GetCentralityPercentile("TRK");
 883     Float_t centv0 = event->GetCentrality()->GetCentralityPercentile("V0M");
 884
 885     if (TMath::Abs(centv0 - centTrk) > 5.0){
 886         centralitypass = kFALSE;
 887         fCentralityNoPass->Fill(fCentrality);
 888      }
 889     const Int_t nGoodTracks = event->GetNumberOfTracks();
 890
 891     Float_t multTPC(0.); // tpc mult estimate
 892     Float_t multGlob(0.); // global multiplicity
 893     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
 894         AliAODTrack* trackAOD = event->GetTrack(iTracks);
 895         if (!trackAOD) continue;
 896         if (!(trackAOD->TestFilterBit(1))) continue;
 897         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;
 898         multTPC++;
 899     }
 900     for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
 901         AliAODTrack* trackAOD = event->GetTrack(iTracks);
 902         if (!trackAOD) continue;
 903         if (!(trackAOD->TestFilterBit(16))) continue;
 904         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;
 905         Double_t b[2] = {-99., -99.};
 906         Double_t bCov[3] = {-99., -99., -99.};
 907         if (!(trackAOD->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
 908         if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
 909         multGlob++;
 910     } //track loop
 911     //     printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
 912  //   if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))){  2010
 913     if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))){
 914         centralitypass = kFALSE;
 915         fCentralityNoPass->Fill(fCentrality);
 916     }//2011
 917     fMultCorBeforeCuts->Fill(multGlob, multTPC);
 918
 919     if(centralitypass){
 920     fCentralityPass->Fill(fCentrality);
 921     fMultCorAfterCuts->Fill(multGlob, multTPC);
 922     fMultvsCentr->Fill(fCentrality, multTPC);
 923     }
 924 }
 925 //_____________________________________________________________________________
 926 void AliAnalysisTaskFlowTPCTOFQCSP::SetCentralityParameters(Double_t CentralityMin, Double_t CentralityMax, const char* CentralityMethod)
 927 {
 928   // Set a centrality range ]min, max] and define the method to use for centrality selection
 929   fCentralityMin = CentralityMin;
 930   fCentralityMax = CentralityMax;
 931   fkCentralityMethod = CentralityMethod;
 932 }
 933 //_____________________________________________________________________________
 934 void AliAnalysisTaskFlowTPCTOFQCSP::SetIDCuts(Double_t minTPCnsigma, Double_t maxTPCnsigma, Double_t minTOFnSigma, Double_t maxTOFnSigma)
 935 {
 936     //Set ID cuts
 937     fminTPCnsigma = minTPCnsigma;
 938     fmaxTPCnsigma = maxTPCnsigma;
 939     fminTOFnSigma = minTOFnSigma;
 940     fmaxTOFnSigma = maxTOFnSigma;
 941 }
 942 //_____________________________________________________________________________
 943
 944 //_____________________________________________________________________________
 945
 946