PWGHF/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2010, 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 // AliAnalysisTaskSEHFv2 gives the needed tools for the D
  19 // mesons v2 analysis with event plane method
  20 // Authors: Chiara Bianchin, cbianchi@pd.infn.it,
  21 // Robert Grajcarek, grajcarek@physi.uni-heidelberg.de
  22 // Giacomo Ortona, ortona@to.infn.it
  23 // Carlos Eugenio Perez Lara, carlos.eugenio.perez.lara@cern.ch
  24 //
  25 /////////////////////////////////////////////////////////////
  26
  27 /* $Id$ */
  28
  29 #include <Riostream.h>
  30 #include <TClonesArray.h>
  31 #include <TCanvas.h>
  32 #include <TList.h>
  33 #include <TFile.h>
  34 #include <TString.h>
  35 #include <TH1F.h>
  36 #include <TH2F.h>
  37 #include <TGraphErrors.h>
  38 #include <TGraph.h>
  39 #include <TDatabasePDG.h>
  40 #include <TRandom3.h>
  41 #include <TVector2.h>
  42 #include <TArrayF.h>
  43
  44 #include <AliLog.h>
  45 #include <AliAnalysisDataSlot.h>
  46 #include <AliAnalysisDataContainer.h>
  47 #include "AliAnalysisManager.h"
  48 #include "AliAODHandler.h"
  49 #include "AliAODEvent.h"
  50 #include "AliAODVertex.h"
  51 #include "AliAODTrack.h"
  52 #include "AliAODMCHeader.h"
  53 #include "AliAODMCParticle.h"
  54 #include "AliAODRecoDecayHF3Prong.h"
  55 #include "AliAODRecoDecayHF.h"
  56 #include "AliAODRecoDecayHF2Prong.h"
  57 #include "AliAODRecoDecayHF4Prong.h"
  58 #include "AliAODRecoCascadeHF.h"
  59
  60 #include "AliAnalysisTaskSE.h"
  61 #include "AliRDHFCutsDplustoKpipi.h"
  62 #include "AliRDHFCutsD0toKpipipi.h"
  63 #include "AliRDHFCutsDstoKKpi.h"
  64 #include "AliRDHFCutsDStartoKpipi.h"
  65 #include "AliRDHFCutsD0toKpi.h"
  66 #include "AliRDHFCutsLctopKpi.h"
  67
  68 #include "AliHFMassFitter.h"
  69 #include "AliEventplane.h"
  70 #include "AliFlowTrack.h"
  71 #include "AliFlowVector.h"
  72 #include "AliFlowTrackCuts.h"
  73 #include "AliFlowEvent.h"
  74
  75 #include "AliAnalysisTaskSEHFv2.h"
  76
  77 ClassImp(AliAnalysisTaskSEHFv2)
  78
  79
  80 //________________________________________________________________________
  81 AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2():
  82 AliAnalysisTaskSE(),
  83   fhEventsInfo(0),
  84   fOutput(0),
  85   fRDCuts(0),
  86   fLowmasslimit(1.669),
  87   fUpmasslimit(2.069),
  88   fNPtBins(1),
  89   fNMassBins(200),
  90   fReadMC(kFALSE),
  91   fUseAfterBurner(kFALSE),
  92   fDecChannel(0),
  93   fAfterBurner(0),
  94   fEventPlaneMeth(kTPCVZERO),
  95   fEventPlanesComp(10),
  96   fV0EPorder(2),
  97   fMinCentr(20),
  98   fMaxCentr(80),
  99   fEtaGap(kFALSE),
 100   fSubEvents(0)
 101 {
 102   // Default constructor
 103 }
 104
 105 //________________________________________________________________________
 106 AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2(const char *name,AliRDHFCuts *rdCuts,Int_t decaychannel):
 107   AliAnalysisTaskSE(name),
 108   fhEventsInfo(0),
 109   fOutput(0),
 110   fRDCuts(rdCuts),
 111   fLowmasslimit(0),
 112   fUpmasslimit(0),
 113   fNPtBins(1),
 114   fNMassBins(200),
 115   fReadMC(kFALSE),
 116   fUseAfterBurner(kFALSE),
 117   fDecChannel(decaychannel),
 118   fAfterBurner(0),
 119   fEventPlaneMeth(kTPCVZERO),
 120   fEventPlanesComp(10),
 121   fV0EPorder(2),
 122   fMinCentr(20),
 123   fMaxCentr(80),
 124   fEtaGap(kFALSE),
 125   fSubEvents(0)
 126 {
 127   // standard constructor
 128   Int_t pdg=421;
 129   switch(fDecChannel){
 130   case 0:
 131     pdg=411;
 132     break;
 133   case 1:
 134     pdg=421;
 135     break;
 136   case 2:
 137     pdg=413;
 138     break;
 139   case 3:
 140     pdg=431;
 141     break;
 142   case 4:
 143     pdg=421;
 144     break;
 145   case 5:
 146     pdg=4122;
 147     break;
 148   }
 149   fAfterBurner = new AliHFAfterBurner(fDecChannel);
 150   if(pdg==413) SetMassLimits((Float_t)0.1,(Float_t)0.2);
 151   else SetMassLimits((Float_t)0.2,pdg); //check range
 152   fNPtBins=fRDCuts->GetNPtBins();
 153
 154   if(fDebug>1)fRDCuts->PrintAll();
 155   // Output slot #1 writes into a TH1F container
 156   DefineOutput(1,TH1F::Class());   //Info on the number of events etc.
 157   // Output slot #2 writes into a TList container
 158   DefineOutput(2,TList::Class());  //Main output
 159   // Output slot #3 writes into a AliRDHFCuts container (cuts)
 160   switch(fDecChannel){
 161   case 0:
 162     DefineOutput(3,AliRDHFCutsDplustoKpipi::Class());  //Cut object for Dplus
 163     break;
 164   case 1:
 165     DefineOutput(3,AliRDHFCutsD0toKpi::Class());  //Cut object for D0
 166     break;
 167   case 2:
 168     DefineOutput(3,AliRDHFCutsDStartoKpipi::Class());  //Cut object for D*
 169     break;
 170   }
 171   //DefineOutput(4,AliFlowEventSimple::Class());
 172   //DefineOutput(4,TList::Class());
 173 }
 174
 175 //________________________________________________________________________
 176 AliAnalysisTaskSEHFv2::~AliAnalysisTaskSEHFv2()
 177 {
 178   // Destructor
 179   delete fOutput;
 180   delete fhEventsInfo;
 181   delete fRDCuts;
 182   delete fAfterBurner;
 183 }
 184 //_________________________________________________________________
 185 void  AliAnalysisTaskSEHFv2::SetMassLimits(Float_t range, Int_t pdg){
 186   // Set limits for mass spectra plots
 187   Float_t mass=0;
 188   Int_t abspdg=TMath::Abs(pdg);
 189   mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
 190   fUpmasslimit = mass+range;
 191   fLowmasslimit = mass-range;
 192 }
 193 //_________________________________________________________________
 194 void  AliAnalysisTaskSEHFv2::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 195   // Set limits for mass spectra plots
 196   if(uplimit>lowlimit)
 197     {
 198       fUpmasslimit = uplimit;
 199       fLowmasslimit = lowlimit;
 200     }
 201 }
 202
 203
 204 //________________________________________________________________________
 205 void AliAnalysisTaskSEHFv2::LocalInit()
 206 {
 207   // Initialization
 208
 209   if(fDebug > 1) printf("AnalysisTaskSEHFv2::Init() \n");
 210
 211   fRDCuts->SetMinCentrality(fMinCentr);
 212   fRDCuts->SetMaxCentrality(fMaxCentr);
 213
 214   switch(fDecChannel){
 215   case 0:
 216     {
 217       AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCuts)));
 218       // Post the data
 219       PostData(3,copycut);
 220     }
 221     break;
 222   case 1:
 223     {
 224       AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCuts)));
 225       // Post the data
 226       PostData(3,copycut);
 227     }
 228     break;
 229   case 2:
 230     {
 231       AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCuts)));
 232       // Post the data
 233       PostData(3,copycut);
 234     }
 235     break;
 236   default:
 237     return;
 238   }
 239   return;
 240 }
 241 //________________________________________________________________________
 242 void AliAnalysisTaskSEHFv2::UserCreateOutputObjects()
 243 {
 244   // Create the output container
 245
 246   if(fDebug > 1) printf("AnalysisTaskSEHFv2::UserCreateOutputObjects() \n");
 247
 248   fhEventsInfo=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Number of AODs scanned",7,-0.5,6.5);
 249   fhEventsInfo->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 250   fhEventsInfo->GetXaxis()->SetBinLabel(2,"nEvSelected");
 251   fhEventsInfo->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
 252   fhEventsInfo->GetXaxis()->SetBinLabel(4,"out of pt bounds");
 253   fhEventsInfo->GetXaxis()->SetBinLabel(5,Form("Ev Sel in Centr %.0f-%.0f%s",fRDCuts->GetMinCentrality(),fRDCuts->GetMaxCentrality(),"%"));
 254   fhEventsInfo->GetXaxis()->SetBinLabel(6,"mismatch lab");
 255   fhEventsInfo->GetXaxis()->SetBinLabel(7,"non valid TPC EP");
 256   fhEventsInfo->GetXaxis()->SetNdivisions(1,kFALSE);
 257
 258
 259   // Several histograms are more conveniently managed in a TList
 260   fOutput = new TList();
 261   fOutput->SetOwner();
 262   fOutput->SetName("MainOutput");
 263
 264   for(Int_t icentr=fMinCentr+5;icentr<=fMaxCentr;icentr=icentr+5){
 265     TString centrname;centrname.Form("centr%d_%d",icentr-5,icentr);
 266
 267       TH2F* hMPtCand=new TH2F(Form("hMPtCand%s",centrname.Data()),Form("Mass vs pt %s;pt (GeV);M (GeV/c^{2})",centrname.Data()),120,0,24.,fNMassBins,fLowmasslimit,fUpmasslimit);
 268       fOutput->Add(hMPtCand);//For <pt> calculation
 269
 270
 271     //Candidate distributions
 272     for(Int_t i=0;i<fNPtBins;i++){
 273       TH2F* hMc2phi=new TH2F(Form("hMc2phi_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
 274       fOutput->Add(hMc2phi);//for 2D analysis
 275
 276       Double_t maxphi = TMath::Pi();
 277       if (fDecChannel == 2) maxphi = TMath::PiOver2();
 278       TH2F* hMphi=new TH2F(Form("hMphi_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi %s;#Delta#phi;M (GeV/c^{2})",centrname.Data()),96,0,maxphi,fNMassBins,fLowmasslimit,fUpmasslimit);
 279       fOutput->Add(hMphi);//for phi bins analysis
 280
 281       if (fReadMC){
 282         TH2F* hMc2phiS=new TH2F(Form("hMc2phiS_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
 283         fOutput->Add(hMc2phiS);
 284         TH2F * hMphiS=new TH2F(Form("hMphiS_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
 285         fOutput->Add(hMphiS);
 286         TH2F* hMc2phiB=new TH2F(Form("hMc2phiB_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
 287         fOutput->Add(hMc2phiB);
 288         TH2F * hMphiB=new TH2F(Form("hMphiB_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
 289         fOutput->Add(hMphiB);
 290         if((fDecChannel != AliAnalysisTaskSEHFv2::kDplustoKpipi) &&(fDecChannel != AliAnalysisTaskSEHFv2::kDstartoKpipi)){
 291           TH2F* hMc2phiR=new TH2F(Form("hMc2phiR_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
 292           fOutput->Add(hMc2phiR);
 293           TH2F* hMphiR=new TH2F(Form("hMphiR_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
 294           fOutput->Add(hMphiR);
 295         }
 296       }
 297     }
 298
 299
 300     //Event Plane
 301     TH2F* hEvPlane=new TH2F(Form("hEvPlane%s",centrname.Data()),Form("VZERO/TPC Event plane angle %s;#phi Ev Plane (TPC);#phi Ev Plane (VZERO);Entries",centrname.Data()),200,0.,TMath::Pi(),200,0.,TMath::Pi());
 302     fOutput->Add(hEvPlane);
 303
 304     TH1F* hEvPlaneA=new TH1F(Form("hEvPlaneA%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
 305     fOutput->Add(hEvPlaneA);
 306
 307     TH1F* hEvPlaneB=new TH1F(Form("hEvPlaneB%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
 308     fOutput->Add(hEvPlaneB);
 309
 310     TH1F* hEvPlaneCand=new TH1F(Form("hEvPlaneCand%s",centrname.Data()),Form("Event plane angle - Event plane angle per candidate %s;#phi(Ev Plane Candidate);Entries",centrname.Data()),200,-TMath::Pi(),TMath::Pi());
 311     fOutput->Add(hEvPlaneCand);
 312
 313     TH1F* hEvPlaneReso=new TH1F(Form("hEvPlaneReso%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{A}-#psi_{B});Entries",centrname.Data()),220,-1.1,1.1);
 314     fOutput->Add(hEvPlaneReso);
 315     if(fEventPlaneMeth>kTPCVZERO){
 316       TH1F* hEvPlaneReso2=new TH1F(Form("hEvPlaneReso2%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{A}-#psi_{B});Entries",centrname.Data()),220,-1.1,1.1);
 317       fOutput->Add(hEvPlaneReso2);
 318       TH1F* hEvPlaneReso3=new TH1F(Form("hEvPlaneReso3%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{A}-#psi_{B});Entries",centrname.Data()),220,-1.1,1.1);
 319       fOutput->Add(hEvPlaneReso3);
 320     }
 321   }
 322
 323   PostData(1,fhEventsInfo);
 324   PostData(2,fOutput);
 325
 326   return;
 327 }
 328
 329 //________________________________________________________________________
 330 void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
 331 {
 332   // Execute analysis for current event:
 333   // heavy flavor candidates association to MC truth
 334   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 335   if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
 336   // Post the data already here
 337   PostData(1,fhEventsInfo);
 338   PostData(2,fOutput);
 339
 340   TClonesArray *arrayProng =0;
 341   Int_t absPdgMom=0;
 342   if(!aod && AODEvent() && IsStandardAOD()) {
 343     // In case there is an AOD handler writing a standard AOD, use the AOD
 344     // event in memory rather than the input (ESD) event.
 345     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 346     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 347     // have to taken from the AOD event hold by the AliAODExtension
 348     AliAODHandler* aodHandler = (AliAODHandler*)
 349       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 350     if(aodHandler->GetExtensions()) {
 351
 352       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 353       AliAODEvent *aodFromExt = ext->GetAOD();
 354
 355
 356       if(fDecChannel==0){
 357         absPdgMom=411;
 358         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 359       }
 360       if(fDecChannel==1){
 361         absPdgMom=421;
 362         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 363       }
 364       if(fDecChannel==2){
 365         absPdgMom=413;
 366         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 367       }
 368     }
 369   } else if(aod){
 370     if(fDecChannel==0){
 371       absPdgMom=411;
 372       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 373     }
 374     if(fDecChannel==1){
 375       absPdgMom=421;
 376       arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 377     }
 378     if(fDecChannel==2){
 379       absPdgMom=413;
 380       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
 381     }
 382   }
 383
 384   if(!aod || !arrayProng) {
 385     AliError("AliAnalysisTaskSEHFv2::UserExec:Branch not found!\n");
 386     return;
 387   }
 388
 389   // fix for temporary bug in ESDfilter
 390   // the AODs with null vertex pointer didn't pass the PhysSel
 391   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 392
 393   TClonesArray *arrayMC=0;
 394   AliAODMCHeader *mcHeader=0;
 395
 396   // load MC particles
 397   if(fReadMC){
 398
 399     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 400     if(!arrayMC) {
 401       AliWarning("AliAnalysisTaskSEHFv2::UserExec:MC particles branch not found!\n");
 402       return;
 403     }
 404
 405     // load MC header
 406     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 407     if(!mcHeader) {
 408       AliError("AliAnalysisTaskSEHFv2::UserExec:MC header branch not found!\n");
 409       return;
 410     }
 411   }
 412
 413   fhEventsInfo->Fill(0); // count event
 414
 415   AliAODRecoDecayHF *d=0;
 416
 417   Int_t nCand = arrayProng->GetEntriesFast();
 418   if(fDebug>2) printf("Number of D2H: %d\n",nCand);
 419
 420   Bool_t isEvSel=fRDCuts->IsEventSelected(aod);
 421   if(!isEvSel){
 422     if(!fRDCuts->IsEventRejectedDueToCentrality())fhEventsInfo->Fill(4);
 423     return;
 424   }
 425
 426   fhEventsInfo->Fill(1);
 427   fhEventsInfo->Fill(4);
 428
 429   AliEventplane *pl=aod->GetEventplane();
 430   if(!pl){
 431     AliError("AliAnalysisTaskSEHFv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
 432     return;
 433   }
 434
 435   //Event plane
 436   Double_t eventplane=0;
 437   Double_t rpangleTPC=0;
 438   Double_t rpangleVZERO=0;
 439   Double_t planereso=0;
 440   Double_t deltaPsi=0;
 441   Double_t rpangleeventB=0;
 442   Double_t rpangleeventA=0;
 443
 444   //For candidate removal from TPC EP
 445   TVector2* q=0x0;
 446   TVector2 *qsub1=0x0;
 447   TVector2 *qsub2=0x0;
 448
 449   //determine centrality bin
 450   Float_t centr=fRDCuts->GetCentrality(aod);
 451   Int_t icentr=0;
 452   for(Int_t ic=fMinCentr+5;ic<=fMaxCentr;ic=ic+5){
 453     if(ic>centr){
 454       icentr=ic;
 455       break;
 456     }
 457   }
 458   TString centrbinname=Form("centr%d_%d",icentr-5,icentr);
 459
 460   if(fReadMC){
 461     TRandom3 *g = new TRandom3(0);
 462     eventplane=g->Rndm()*TMath::Pi();
 463     delete g;g=0x0;
 464     if(fUseAfterBurner)fAfterBurner->SetEventPlane((Double_t)eventplane);
 465   }else{
 466     if(fEventPlaneMeth!=kTPC){
 467
 468       //VZERO EP and resolution
 469       rpangleVZERO=GetPhi0Pi(pl->GetEventplane("V0",aod,fV0EPorder));
 470       if(fEventPlaneMeth>kTPCVZERO){
 471         //Using V0A/C for VZERO resolution
 472         if(fEventPlaneMeth==kVZEROA){
 473           rpangleVZERO=GetPhi0Pi(pl->GetEventplane("V0A",aod,fV0EPorder));
 474           if(fSubEvents!=kSingleV0Side)rpangleeventB=GetPhi0Pi(pl->GetEventplane("V0C",aod,fV0EPorder));
 475         }
 476         else if(fEventPlaneMeth==kVZEROC){
 477           rpangleVZERO=GetPhi0Pi(pl->GetEventplane("V0C",aod,fV0EPorder));
 478           if(fSubEvents!=kSingleV0Side)rpangleeventB=GetPhi0Pi(pl->GetEventplane("V0A",aod,fV0EPorder));
 479         }
 480         // deltaPsi =rpangleeventA-rpangleeventB;
 481         eventplane=rpangleVZERO;
 482       }
 483     }
 484
 485     // TPC event plane
 486     rpangleTPC = pl->GetEventplane("Q");
 487     if(rpangleTPC<0){
 488       fhEventsInfo->Fill(6);
 489       return;
 490     }
 491     rpangleeventA = rpangleTPC;
 492     if(fSubEvents>kFullTPC||fEtaGap||fEventPlaneMeth==kVZERO){
 493       qsub1 = pl->GetQsub1();
 494       qsub2 = pl->GetQsub2();
 495       if(!qsub1 || !qsub2){
 496         fhEventsInfo->Fill(6);
 497         return;
 498       }
 499       if(fSubEvents==kPosTPC||fSubEvents==kSingleV0Side||fEventPlaneMeth==kVZERO){
 500         rpangleeventA = qsub1->Phi()/2.;
 501         if(fSubEvents==kSingleV0Side||fEventPlaneMeth==kVZERO)rpangleeventB = qsub2->Phi()/2.;
 502       }
 503       if(fSubEvents==kNegTPC)rpangleeventA = qsub2->Phi()/2.;
 504     }
 505
 506     if(fEventPlaneMeth<=kTPCVZERO){
 507       q = pl->GetQVector();
 508       deltaPsi = pl->GetQsubRes();
 509     }else{
 510       deltaPsi=eventplane-rpangleeventA;
 511     }
 512   }
 513
 514   if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){
 515     if(deltaPsi>0.) deltaPsi-=TMath::Pi();
 516     else deltaPsi +=TMath::Pi();
 517   } // difference of subevents reaction plane angle cannot be bigger than phi/2
 518   planereso = TMath::Cos(2.*deltaPsi); // reaction plane resolution
 519   if(TMath::Abs(rpangleTPC-rpangleVZERO)>fEventPlanesComp)return;
 520
 521   if(fDebug>2)printf("Filling EP-related histograms\n");
 522   //Filling EP-related histograms
 523   ((TH2F*)fOutput->FindObject(Form("hEvPlane%s",centrbinname.Data())))->Fill(rpangleTPC,rpangleVZERO); // reaction plane angle without autocorrelations removal
 524   ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso%s",centrbinname.Data())))->Fill(planereso); //RP resolution
 525   if(fEventPlaneMeth>kTPCVZERO||fEtaGap){
 526     ((TH1F*)fOutput->FindObject(Form("hEvPlaneA%s",centrbinname.Data())))->Fill(rpangleeventA); //Angle of first subevent
 527     ((TH1F*)fOutput->FindObject(Form("hEvPlaneB%s",centrbinname.Data())))->Fill(rpangleeventB); //Angle of second subevent
 528   }
 529   if(fEventPlaneMeth>kTPCVZERO){
 530     Double_t deltaSub=rpangleeventA-rpangleeventB;
 531     if(TMath::Abs(deltaSub)>TMath::Pi()/2.){// difference of subevents reaction plane angle cannot be bigger than phi/2
 532       if(deltaSub>0.) deltaSub-=TMath::Pi();
 533       else deltaSub +=TMath::Pi();
 534     }
 535     ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso2%s",centrbinname.Data())))->Fill(TMath::Cos(2.*deltaSub)); //RP resolution
 536     deltaSub =eventplane-rpangleeventB;
 537     if(TMath::Abs(deltaSub)>TMath::Pi()/2.){// difference of subevents reaction plane angle cannot be bigger than phi/2
 538       if(deltaSub>0.) deltaSub-=TMath::Pi();
 539       else deltaSub +=TMath::Pi();
 540     }
 541     ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso3%s",centrbinname.Data())))->Fill(TMath::Cos(2.*deltaSub)); //RP resolution
 542   }
 543
 544   if(fDebug>2)printf("Loop on D candidates\n");
 545   //Loop on D candidates
 546   for (Int_t iCand = 0; iCand < nCand; iCand++) {
 547     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iCand);
 548     Bool_t isSelBit=kTRUE;
 549     if(fDecChannel==0) isSelBit=d->HasSelectionBit(AliRDHFCuts::kDplusCuts);
 550     if(fDecChannel==1) isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts);
 551     if(fDecChannel==2) isSelBit=kTRUE;
 552     if(!isSelBit)continue;
 553     Int_t ptbin=fRDCuts->PtBin(d->Pt());
 554     if(ptbin<0) {
 555       fhEventsInfo->Fill(3);
 556       continue;
 557     }
 558     Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(absPdgMom));
 559     if(!isFidAcc)continue;
 560     Int_t isSelected= fRDCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 561     if(!isSelected)continue;
 562
 563     fhEventsInfo->Fill(2); // candidate selected
 564     if(fDebug>3) printf("+++++++Is Selected\n");
 565
 566     Float_t* invMass=0x0;
 567     Int_t nmasses;
 568     CalculateInvMasses(d,invMass,nmasses);
 569
 570     if(fEventPlaneMeth<=kTPCVZERO){
 571       eventplane = GetEventPlaneForCandidate(d,q,pl,qsub1,qsub2); // remove autocorrelations
 572       ((TH1F*)fOutput->FindObject(Form("hEvPlaneCand%s",centrbinname.Data())))->Fill(rpangleTPC-eventplane);
 573     }
 574
 575     Double_t phi=d->Phi();
 576     if(fReadMC&&fUseAfterBurner)phi=fAfterBurner->GetNewAngle(d,arrayMC);
 577     Float_t deltaphi=GetPhi0Pi(phi-eventplane);
 578
 579     //fill the histograms with the appropriate method
 580     if(fDecChannel==0)FillDplus(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
 581     else if(fDecChannel==1)FillD02p(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
 582     else if(fDecChannel==2)FillDstar(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
 583
 584     delete [] invMass;
 585   }
 586
 587   PostData(1,fhEventsInfo);
 588   PostData(2,fOutput);
 589
 590   return;
 591 }
 592
 593 //***************************************************************************
 594
 595 // Methods used in the UserExec
 596
 597 void AliAnalysisTaskSEHFv2::CalculateInvMasses(AliAODRecoDecayHF* d,Float_t*& masses,Int_t& nmasses){
 598   //Calculates all the possible invariant masses for each candidate
 599   //NB: the implementation for each candidate is responsibility of the corresponding developer
 600
 601   if(fDecChannel==0){
 602     //D+ -- Giacomo
 603     nmasses=1;
 604     masses=new Float_t[nmasses];
 605     Int_t pdgdaughters[3] = {211,321,211};
 606     masses[0]=d->InvMass(3,(UInt_t*)pdgdaughters);
 607   }
 608   if(fDecChannel==1){
 609     //D0 (Kpi)  -- Chiara
 610     const Int_t ndght=2;
 611     nmasses=2;
 612     masses=new Float_t[nmasses];
 613     Int_t pdgdaughtersD0[ndght]={211,321};//pi,K
 614     masses[0]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0); //D0
 615     Int_t pdgdaughtersD0bar[ndght]={321,211};//K,pi
 616     masses[1]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0bar); //D0bar
 617   }
 618   if(fDecChannel==2){
 619     //D* -- Robert,Yifei, Alessandro
 620     nmasses=1;
 621     masses=new Float_t[nmasses];
 622     masses[0]=((AliAODRecoCascadeHF*)d)->DeltaInvMass();
 623   }
 624 }
 625
 626 //******************************************************************************
 627
 628 //Methods to fill the histograms, one for each channel
 629 //NB: the implementation for each candidate is responsibility of the corresponding developer
 630
 631 //******************************************************************************
 632 void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
 633   //D+ channel
 634   if(!isSel){
 635     if(fDebug>3)AliWarning("Candidate not selected\n");
 636     return;
 637   }
 638   if(!masses){
 639     if(fDebug>3)AliWarning("Masses not calculated\n");
 640     return;
 641   }
 642
 643   ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 644   ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 645   ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
 646   Int_t pdgdaughters[3] = {211,321,211};
 647
 648   if(fReadMC){
 649     Int_t lab=-1;
 650     if(fUseAfterBurner){
 651       Bool_t isSignal=fAfterBurner->GetIsSignal();
 652       if(isSignal)lab=10;
 653     }else {
 654       lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
 655     }
 656     if(lab>=0){ //signal
 657       ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 658       ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 659     } else{ //background
 660       ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 661       ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 662     }
 663   }
 664 }
 665
 666 //******************************************************************************
 667 void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
 668
 669   //D0->Kpi channel
 670
 671   //mass histograms
 672   if(!masses){
 673     if(fDebug>3)AliWarning("Masses not calculated\n");
 674     return;
 675   }
 676   if(isSel==1 || isSel==3) {
 677     ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 678     ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 679     ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
 680   }
 681   if(isSel>=2) {
 682     ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 683     ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
 684     ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[1]);
 685   }
 686
 687   //MC histograms
 688   if(fReadMC){
 689
 690     Int_t matchtoMC=-1;
 691
 692     //D0
 693     Int_t pdgdaughters[2];
 694     pdgdaughters[0]=211;//pi
 695     pdgdaughters[1]=321;//K
 696     Int_t nprongs=2;
 697     Int_t absPdgMom=421;
 698
 699     matchtoMC = d->MatchToMC(absPdgMom,arrayMC,nprongs,pdgdaughters);
 700
 701     Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
 702     if((isSel==1 || isSel==3)){ //D0
 703       if(matchtoMC>=0){
 704         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 705         Int_t pdgMC = dMC->GetPdgCode();
 706
 707         if(pdgMC==prongPdgPlus) {
 708           ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 709           ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 710         }
 711         else {
 712           ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 713           ((TH2F*)fOutput->FindObject(Form("hMphiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 714         }
 715       } else {
 716         ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 717         ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 718       }
 719     }
 720     if(isSel>=2){ //D0bar
 721       if(matchtoMC>=0){
 722         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 723         Int_t pdgMC = dMC->GetPdgCode();
 724
 725         if(pdgMC==prongPdgMinus) {
 726           ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 727           ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
 728         }
 729         else {
 730           ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 731           ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
 732         }
 733       } else {
 734         ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 735         ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
 736       }
 737     }
 738   }
 739 }
 740 //******************************************************************************
 741 void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
 742   //D* channel
 743   if(!isSel){
 744     if(fDebug>3)AliWarning("Candidate not selected\n");
 745     return;
 746   }
 747   if(!masses){
 748     if(fDebug>3)AliWarning("Masses not calculated\n");
 749     return;
 750   }
 751
 752   Float_t deltaphi1 = deltaphi;
 753   if(deltaphi1 > TMath::PiOver2()) deltaphi1 = TMath::Pi() - deltaphi1;
 754
 755   ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 756   ((TH2F*)fOutput->FindObject(Form("hMphi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi1,masses[0]);
 757   ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
 758   Int_t pdgDgDStartoD0pi[2]={421,211};
 759   Int_t pdgDgD0toKpi[2]={321,211};
 760
 761   if(fReadMC){
 762     Int_t lab=-1;
 763     lab = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,arrayMC);
 764     if(lab>=0){ //signal
 765       ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 766       ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 767     } else{ //background
 768       ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 769       ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 770     }
 771   }
 772 }
 773
 774 //________________________________________________________________________
 775 void AliAnalysisTaskSEHFv2::SetEventPlaneMethod(Int_t method){
 776   if(method>kVZEROC||method<0){
 777     AliWarning("No EP method associated to the selection, setting to TPC EP\n");
 778     method=kTPCVZERO;
 779   }
 780   fEventPlaneMeth=method;
 781 }
 782
 783 //________________________________________________________________________
 784 Float_t AliAnalysisTaskSEHFv2::GetPhi0Pi(Float_t phi){
 785   // Sets the phi angle in the range 0-pi
 786   Float_t result=phi;
 787   while(result<0){
 788     result=result+TMath::Pi();
 789   }
 790   while(result>TMath::Pi()){
 791     result=result-TMath::Pi();
 792   }
 793    return result;
 794 }
 795
 796 //________________________________________________________________________
 797 Float_t AliAnalysisTaskSEHFv2::GetEventPlaneForCandidate(AliAODRecoDecayHF* d, const TVector2* q,AliEventplane *pl, const TVector2* qsub1, const TVector2* qsub2){
 798   // remove autocorrelations
 799
 800   TArrayF* qx = 0x0;
 801   TArrayF* qy = 0x0;
 802   TVector2 qcopy;
 803   if(!fEtaGap){
 804     qx = pl->GetQContributionXArray();
 805     qy = pl->GetQContributionYArray();
 806     qcopy = *q;
 807     }
 808   else {
 809     if(d->Eta()>0.){
 810       qx = pl->GetQContributionXArraysub1();
 811       qy = pl->GetQContributionYArraysub1();
 812       qcopy = *qsub1;
 813     }
 814     else{
 815       qx = pl->GetQContributionXArraysub2();
 816       qy = pl->GetQContributionYArraysub2();
 817       qcopy = *qsub2;
 818     }
 819   }
 820
 821   if(fDecChannel==2){
 822     //D* -- Yifei, Alessandro,Robert
 823     AliAODRecoDecayHF2Prong* theD0particle = ((AliAODRecoCascadeHF*)d)->Get2Prong();
 824     AliAODTrack *track0 = (AliAODTrack*)theD0particle->GetDaughter(0);
 825     AliAODTrack *track1 = (AliAODTrack*)theD0particle->GetDaughter(1);
 826     AliAODTrack *track2 = ((AliAODRecoCascadeHF*)d)->GetBachelor();
 827     // reduce global q vector
 828
 829     TVector2 q0;
 830     if((track0->GetID()) < qx->fN){
 831       q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
 832
 833     TVector2 q1;
 834     if((track1->GetID()) < qx->fN){
 835       q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
 836
 837     TVector2 q2;
 838     if((track2->GetID()) < qx->fN){
 839       q2.Set(qx->At(track2->GetID()),qy->At(track2->GetID()));}
 840
 841     qcopy = qcopy -(q0+q1+q2);
 842
 843   }
 844
 845   // reduce Q vector for D+ and D0
 846
 847   if(fDecChannel==1){
 848     //D0 -- Chiara
 849     AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
 850     AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
 851
 852     TVector2 q0;
 853     if((track0->GetID()) < qx->fN){
 854       q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
 855
 856     TVector2 q1;
 857     if((track1->GetID()) < qx->fN){
 858       q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
 859
 860     qcopy = qcopy -(q0+q1);
 861   }
 862
 863   if(fDecChannel==0){
 864     //D+ -- Giacomo
 865     AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
 866     AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
 867     AliAODTrack *track2 = (AliAODTrack*)d->GetDaughter(2);
 868
 869     TVector2 q0;
 870     if((track0->GetID()) < qx->fN){
 871       q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
 872
 873     TVector2 q1;
 874     if((track1->GetID()) < qx->fN){
 875       q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
 876
 877     TVector2 q2;
 878     if((track2->GetID()) < qx->fN){
 879       q2.Set(qx->At(track2->GetID()),qy->At(track2->GetID()));}
 880
 881     qcopy = qcopy -(q0+q1+q2);
 882
 883   }
 884
 885   return qcopy.Phi()/2.;
 886
 887 }
 888 // //________________________________________________________________________
 889 // Float_t AliAnalysisTaskSEHFv2::GetEventPlaneFromV0(AliAODEvent *aodEvent){
 890 //   // Compute event plane for VZERO - Obsolete, used for 2010 data
 891
 892 //   Int_t centr=fRDCuts->GetCentrality(aodEvent);
 893 //   centr=centr-centr%10;
 894 //   //temporary fix
 895 //   if(centr<20)centr=20;
 896 //   if(centr>70)centr=70;
 897 //   //end temporary fix
 898 //   Int_t binx=0;
 899 //   Int_t iParHist=(centr-20)/10;
 900
 901 //   TString name;name.Form("parhist%d_%d",centr,centr+10);
 902
 903 //   if(fDebug>15)printf("EPfromV0 centr %d, iparhist %d (%p-%p)\n",centr,iParHist,fParHist->FindObject(name.Data()),fParHist->At(iParHist));
 904
 905 //   Int_t runnumber=aodEvent->GetRunNumber();
 906 //   if(fParHist->At(iParHist)){
 907 //     for(Int_t i=1;i<=((TH2D*)fParHist->At(iParHist))->GetNbinsX()&&binx<=0;i++){
 908 //       Int_t run=atoi(((TH2D*)fParHist->At(iParHist))->GetXaxis()->GetBinLabel(i));
 909 //       if(run>=runnumber)binx=i;
 910 //     }
 911 //   }else{
 912 //     fhEventsInfo->Fill(7);
 913 //   }
 914
 915 //   AliFlowTrackCuts* cutsRP = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts();
 916 //   cutsRP->SetEvent(aodEvent, MCEvent());//, 0x0);
 917 //   cutsRP->SetName( Form("rp_cuts") );
 918 //   AliFlowTrackCuts* dummy = new AliFlowTrackCuts("null_cuts");
 919 //   dummy->SetParamType(AliFlowTrackCuts::kGlobal);
 920 //   dummy->SetPtRange(+1,-1); // select nothing QUICK
 921 //   dummy->SetEtaRange(+1,-1); // select nothing VZERO
 922 //   dummy->SetEvent(aodEvent,MCEvent());
 923
 924 //   //////////////// construct the flow event container ////////////
 925 //   AliFlowEvent flowEvent(cutsRP,dummy);
 926 //   flowEvent.SetReferenceMultiplicity( 64 );
 927 //   for(Int_t i=0;i<64&&binx>0;i++){
 928 //     AliFlowTrack *flowTrack=flowEvent.GetTrack(i);
 929 //     Double_t inte=((TH2D*)fParHist->At(iParHist))->Integral(binx,binx,i+1,i+1);
 930 //     if(inte>0)flowTrack->SetWeight(flowTrack->Weight()/inte);
 931 //   }
 932 //   if(fDebug>15)printf("EPfromV0 flow tracks weights done\n");
 933
 934 //   AliFlowVector qvec=flowEvent.GetQ(fV0EPorder);
 935 //   Double_t angleEP=(1./(Double_t)fV0EPorder)*qvec.Phi();
 936 //   if(fDebug>15)printf("EPfromV0 phi %f\n",angleEP);
 937 //   return angleEP;
 938 // }
 939 //________________________________________________________________________
 940 void AliAnalysisTaskSEHFv2::Terminate(Option_t */*option*/)
 941 {
 942   // Terminate analysis
 943   //
 944   if(fDebug > 1) printf("AnalysisTaskSEHFv2: Terminate() \n");
 945
 946   return;
 947 }
 948