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