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
 275     for(Int_t i=0;i<fNPtBins;i++){
 276
 277       // Delta Phi histograms
 278       Double_t maxphi = TMath::Pi();
 279       if (fDecChannel == 2) maxphi = TMath::PiOver2();
 280       TH2F* hMdeltaphi=new TH2F(Form("hMdeltaphi_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(hMdeltaphi);//for phi bins analysis
 282       TH2F* hMc2deltaphi=new TH2F(Form("hMc2deltaphi_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(hMc2deltaphi);
 284       TH2F* hMs2deltaphi=new TH2F(Form("hMs2deltaphi_pt%d%s",i,centrname.Data()),Form("Mass vs sin2#Delta#phi (p_{t} bin %d %s);sin2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
 285       fOutput->Add(hMs2deltaphi);
 286
 287       // phi histograms (for systematics)
 288       TH2F* hCos2PhiMass=new TH2F(Form("hCos2phiMass_pt%d%s",i,centrname.Data()),Form("Mass vs cos(2#phi) %s;cos(2#phi);M (GeV/c^{2})",centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
 289       fOutput->Add(hCos2PhiMass);
 290       TH2F* hSin2PhiMass=new TH2F(Form("hSin2phiMass_pt%d%s",i,centrname.Data()),Form("Mass vs sin(2#phi) %s;sin(2#phi);M (GeV/c^{2})",centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
 291       fOutput->Add(hSin2PhiMass);
 292
 293       // Histos using MC truth
 294       if (fReadMC){
 295         TH2F* hMc2deltaphiS=new TH2F(Form("hMc2deltaphiS_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);
 296         fOutput->Add(hMc2deltaphiS);
 297         TH2F * hMdeltaphiS=new TH2F(Form("hMdeltaphiS_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);
 298         fOutput->Add(hMdeltaphiS);
 299         TH2F* hMc2deltaphiB=new TH2F(Form("hMc2deltaphiB_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);
 300         fOutput->Add(hMc2deltaphiB);
 301         TH2F * hMdeltaphiB=new TH2F(Form("hMdeltaphiB_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);
 302         fOutput->Add(hMdeltaphiB);
 303         if((fDecChannel != AliAnalysisTaskSEHFv2::kDplustoKpipi) &&(fDecChannel != AliAnalysisTaskSEHFv2::kDstartoKpipi)){
 304           TH2F* hMc2deltaphiR=new TH2F(Form("hMc2deltaphiR_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);
 305           fOutput->Add(hMc2deltaphiR);
 306           TH2F* hMdeltaphiR=new TH2F(Form("hMdeltaphiR_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);
 307           fOutput->Add(hMdeltaphiR);
 308         }
 309       }
 310     }
 311
 312
 313     // Event Plane
 314     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());
 315     fOutput->Add(hEvPlane);
 316     TH1F* hEvPlaneA=new TH1F(Form("hEvPlaneA%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
 317     fOutput->Add(hEvPlaneA);
 318     TH1F* hEvPlaneB=new TH1F(Form("hEvPlaneB%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
 319     fOutput->Add(hEvPlaneB);
 320     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());
 321     fOutput->Add(hEvPlaneCand);
 322
 323     // histos for EP resolution
 324     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);
 325     fOutput->Add(hEvPlaneReso);
 326     if(fEventPlaneMeth>=kTPCVZERO){
 327       TH1F* hEvPlaneReso2=new TH1F(Form("hEvPlaneReso2%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{A}-#psi_{C});Entries",centrname.Data()),220,-1.1,1.1);
 328       fOutput->Add(hEvPlaneReso2);
 329       TH1F* hEvPlaneReso3=new TH1F(Form("hEvPlaneReso3%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{B}-#psi_{C});Entries",centrname.Data()),220,-1.1,1.1);
 330       fOutput->Add(hEvPlaneReso3);
 331     }
 332
 333     // histos for EPsystematics
 334     TH1F *hCos2EP=new TH1F(Form("hCos2EP%s",centrname.Data()),Form("cos(2PsiEP) %s;cos2(#psi_{EP});Entries",centrname.Data()),100,-1.,1.);
 335     fOutput->Add(hCos2EP);
 336     TH1F *hSin2EP=new TH1F(Form("hSin2EP%s",centrname.Data()),Form("sin(2PsiEP) %s;sin2(#psi_{EP});Entries",centrname.Data()),100,-1.,1.);
 337     fOutput->Add(hSin2EP);
 338    }
 339
 340   PostData(1,fhEventsInfo);
 341   PostData(2,fOutput);
 342
 343   return;
 344 }
 345
 346 //________________________________________________________________________
 347 void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
 348 {
 349   // Execute analysis for current event:
 350   // heavy flavor candidates association to MC truth
 351   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 352   if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
 353   // Post the data already here
 354   PostData(1,fhEventsInfo);
 355   PostData(2,fOutput);
 356
 357   TClonesArray *arrayProng =0;
 358   Int_t absPdgMom=0;
 359   if(!aod && AODEvent() && IsStandardAOD()) {
 360     // In case there is an AOD handler writing a standard AOD, use the AOD
 361     // event in memory rather than the input (ESD) event.
 362     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 363     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 364     // have to taken from the AOD event hold by the AliAODExtension
 365     AliAODHandler* aodHandler = (AliAODHandler*)
 366       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 367     if(aodHandler->GetExtensions()) {
 368
 369       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 370       AliAODEvent *aodFromExt = ext->GetAOD();
 371
 372
 373       if(fDecChannel==0){
 374         absPdgMom=411;
 375         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 376       }
 377       if(fDecChannel==1){
 378         absPdgMom=421;
 379         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 380       }
 381       if(fDecChannel==2){
 382         absPdgMom=413;
 383         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 384       }
 385     }
 386   } else if(aod){
 387     if(fDecChannel==0){
 388       absPdgMom=411;
 389       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 390     }
 391     if(fDecChannel==1){
 392       absPdgMom=421;
 393       arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 394     }
 395     if(fDecChannel==2){
 396       absPdgMom=413;
 397       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
 398     }
 399   }
 400
 401   if(!aod || !arrayProng) {
 402     AliError("AliAnalysisTaskSEHFv2::UserExec:Branch not found!\n");
 403     return;
 404   }
 405
 406   // fix for temporary bug in ESDfilter
 407   // the AODs with null vertex pointer didn't pass the PhysSel
 408   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 409
 410   TClonesArray *arrayMC=0;
 411   AliAODMCHeader *mcHeader=0;
 412
 413   // load MC particles
 414   if(fReadMC){
 415
 416     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 417     if(!arrayMC) {
 418       AliWarning("AliAnalysisTaskSEHFv2::UserExec:MC particles branch not found!\n");
 419       return;
 420     }
 421
 422     // load MC header
 423     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 424     if(!mcHeader) {
 425       AliError("AliAnalysisTaskSEHFv2::UserExec:MC header branch not found!\n");
 426       return;
 427     }
 428   }
 429
 430   fhEventsInfo->Fill(0); // count event
 431
 432   AliAODRecoDecayHF *d=0;
 433
 434   Int_t nCand = arrayProng->GetEntriesFast();
 435   if(fDebug>2) printf("Number of D2H: %d\n",nCand);
 436
 437   Bool_t isEvSel=fRDCuts->IsEventSelected(aod);
 438   if(!isEvSel){
 439     if(!fRDCuts->IsEventRejectedDueToCentrality())fhEventsInfo->Fill(4);
 440     return;
 441   }
 442
 443   fhEventsInfo->Fill(1);
 444   fhEventsInfo->Fill(4);
 445
 446   AliEventplane *pl=aod->GetEventplane();
 447   if(!pl){
 448     AliError("AliAnalysisTaskSEHFv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
 449     fhEventsInfo->Fill(6);
 450     return;
 451   }
 452
 453   //Event plane
 454   Double_t eventplane=0;
 455   Double_t rpangleTPC=0;
 456   Double_t rpangleVZERO=0;
 457   Double_t planereso=0;
 458   Double_t deltaPsi=0;
 459   Double_t rpangleeventC=0;
 460   Double_t rpangleeventB=0;
 461   Double_t rpangleeventA=0;
 462
 463   //For candidate removal from TPC EP
 464   TVector2* q=0x0;
 465   TVector2 *qsub1=0x0;
 466   TVector2 *qsub2=0x0;
 467
 468   //determine centrality bin
 469   Float_t centr=fRDCuts->GetCentrality(aod);
 470   Float_t centrPerMil=centr*10.;
 471   Int_t icentr=-1;
 472   for(Int_t ic=fMinCentr*10+fCentBinSizePerMil;ic<=fMaxCentr*10;ic=ic+fCentBinSizePerMil){
 473     if(ic>centrPerMil){
 474       icentr=ic;
 475       break;
 476     }
 477   }
 478   if(icentr==-1) return;
 479
 480   TString centrbinname=Form("centr%d_%d",icentr-fCentBinSizePerMil,icentr);
 481
 482   if(fReadMC){
 483     TRandom3 *g = new TRandom3(0);
 484     eventplane=g->Rndm()*TMath::Pi();
 485     delete g;g=0x0;
 486     if(fUseAfterBurner)fAfterBurner->SetEventPlane((Double_t)eventplane);
 487   }else{
 488    // TPC event plane
 489     rpangleTPC = pl->GetEventplane("Q");
 490     if(rpangleTPC<0){
 491       fhEventsInfo->Fill(6);
 492       return;
 493     }
 494     rpangleeventA = rpangleTPC;
 495     if(fSubEvents==2||fEventPlaneMeth==kVZERO){
 496       qsub1 = pl->GetQsub1();
 497       qsub2 = pl->GetQsub2();
 498       if(!qsub1 || !qsub2){
 499         fhEventsInfo->Fill(6);
 500         return;
 501       }
 502       rpangleeventA = qsub1->Phi()/2.;
 503       rpangleeventB = qsub2->Phi()/2.;
 504     }
 505     if(fEventPlaneMeth!=kTPC){
 506       //VZERO EP and resolution
 507       rpangleVZERO=GetPhi0Pi(pl->GetEventplane("V0",aod,fV0EPorder));
 508       rpangleeventC=rpangleVZERO;
 509       if(fEventPlaneMeth==kVZEROA||fEventPlaneMeth==kVZEROC||(fEventPlaneMeth==kTPCVZERO&&fSubEvents==3)){
 510         rpangleeventB=GetPhi0Pi(pl->GetEventplane("V0A",aod,fV0EPorder));
 511         rpangleeventC=GetPhi0Pi(pl->GetEventplane("V0C",aod,fV0EPorder));
 512         if(fEventPlaneMeth==kVZEROA)rpangleVZERO=rpangleeventB;
 513         else if(fEventPlaneMeth==kVZEROC)rpangleVZERO=rpangleeventC;
 514       }
 515     }
 516
 517
 518     if(fEventPlaneMeth>kTPCVZERO)eventplane=rpangleVZERO;
 519     else{
 520       q = pl->GetQVector();
 521       eventplane=rpangleTPC;
 522     }
 523     deltaPsi =rpangleeventA-rpangleeventB;
 524   }
 525
 526   if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){
 527     if(deltaPsi>0.) deltaPsi-=TMath::Pi();
 528     else deltaPsi +=TMath::Pi();
 529   } // difference of subevents reaction plane angle cannot be bigger than phi/2
 530   planereso = TMath::Cos(2.*deltaPsi); // reaction plane resolution
 531   if(TMath::Abs(rpangleTPC-rpangleVZERO)>fEventPlanesComp)return;
 532
 533   if(fDebug>2)printf("Filling EP-related histograms\n");
 534   //Filling EP-related histograms
 535   ((TH2F*)fOutput->FindObject(Form("hEvPlane%s",centrbinname.Data())))->Fill(rpangleTPC,rpangleVZERO); // reaction plane angle without autocorrelations removal
 536   ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso%s",centrbinname.Data())))->Fill(planereso); //RP resolution
 537   ((TH1F*)fOutput->FindObject(Form("hCos2EP%s",centrbinname.Data())))->Fill(TMath::Cos(2*eventplane)); // syst check
 538   ((TH1F*)fOutput->FindObject(Form("hSin2EP%s",centrbinname.Data())))->Fill(TMath::Sin(2*eventplane)); // syst check
 539
 540   if(fEventPlaneMeth>kTPCVZERO||fEtaGap){
 541     ((TH1F*)fOutput->FindObject(Form("hEvPlaneA%s",centrbinname.Data())))->Fill(rpangleeventA); //Angle of first subevent
 542     ((TH1F*)fOutput->FindObject(Form("hEvPlaneB%s",centrbinname.Data())))->Fill(rpangleeventB); //Angle of second subevent
 543   }
 544   if(fEventPlaneMeth>kTPCVZERO||fSubEvents==3){
 545     Double_t deltaSub=rpangleeventA-rpangleeventC;
 546     if(TMath::Abs(deltaSub)>TMath::Pi()/2.){// difference of subevents reaction plane angle cannot be bigger than phi/2
 547       if(deltaSub>0.) deltaSub-=TMath::Pi();
 548       else deltaSub +=TMath::Pi();
 549     }
 550     TH1F* htmp1=(TH1F*)fOutput->FindObject(Form("hEvPlaneReso2%s",centrbinname.Data()));
 551     if(htmp1) htmp1->Fill(TMath::Cos(2.*deltaSub)); //RP resolution
 552     deltaSub =rpangleeventB-rpangleeventC;
 553     if(TMath::Abs(deltaSub)>TMath::Pi()/2.){// difference of subevents reaction plane angle cannot be bigger than phi/2
 554       if(deltaSub>0.) deltaSub-=TMath::Pi();
 555       else deltaSub +=TMath::Pi();
 556     }
 557     TH1F* htmp2=(TH1F*)fOutput->FindObject(Form("hEvPlaneReso3%s",centrbinname.Data()));
 558     if(htmp2) htmp2->Fill(TMath::Cos(2.*deltaSub)); //RP resolution
 559   }
 560
 561   if(fDebug>2)printf("Loop on D candidates\n");
 562   //Loop on D candidates
 563   for (Int_t iCand = 0; iCand < nCand; iCand++) {
 564     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iCand);
 565     Bool_t isSelBit=kTRUE;
 566     if(fDecChannel==0) isSelBit=d->HasSelectionBit(AliRDHFCuts::kDplusCuts);
 567     if(fDecChannel==1) isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts);
 568     if(fDecChannel==2) isSelBit=kTRUE;
 569     if(!isSelBit)continue;
 570     Int_t ptbin=fRDCuts->PtBin(d->Pt());
 571     if(ptbin<0) {
 572       fhEventsInfo->Fill(3);
 573       continue;
 574     }
 575     Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(absPdgMom));
 576     if(!isFidAcc)continue;
 577     Int_t isSelected= fRDCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
 578     if(!isSelected)continue;
 579
 580     fhEventsInfo->Fill(2); // candidate selected
 581     if(fDebug>3) printf("+++++++Is Selected\n");
 582
 583     Float_t* invMass=0x0;
 584     Int_t nmasses;
 585     CalculateInvMasses(d,invMass,nmasses);
 586
 587     if(fEventPlaneMeth<=kTPCVZERO){
 588       eventplane = GetEventPlaneForCandidate(d,q,pl,qsub1,qsub2); // remove autocorrelations
 589       ((TH1F*)fOutput->FindObject(Form("hEvPlaneCand%s",centrbinname.Data())))->Fill(rpangleTPC-eventplane);
 590     }
 591
 592     Double_t phi=d->Phi();
 593     if(fReadMC&&fUseAfterBurner)phi=fAfterBurner->GetNewAngle(d,arrayMC);
 594     Float_t deltaphi=GetPhi0Pi(phi-eventplane);
 595
 596     //fill the histograms with the appropriate method
 597     if(fDecChannel==0)FillDplus(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr,phi);
 598     else if(fDecChannel==1)FillD02p(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr,phi);
 599     else if(fDecChannel==2)FillDstar(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr,phi);
 600
 601     delete [] invMass;
 602   }
 603
 604   PostData(1,fhEventsInfo);
 605   PostData(2,fOutput);
 606
 607   return;
 608 }
 609
 610 //***************************************************************************
 611
 612 // Methods used in the UserExec
 613
 614 void AliAnalysisTaskSEHFv2::CalculateInvMasses(AliAODRecoDecayHF* d,Float_t*& masses,Int_t& nmasses){
 615   //Calculates all the possible invariant masses for each candidate
 616   //NB: the implementation for each candidate is responsibility of the corresponding developer
 617
 618   if(fDecChannel==0){
 619     //D+ -- Giacomo
 620     nmasses=1;
 621     masses=new Float_t[nmasses];
 622     Int_t pdgdaughters[3] = {211,321,211};
 623     masses[0]=d->InvMass(3,(UInt_t*)pdgdaughters);
 624   }
 625   if(fDecChannel==1){
 626     //D0 (Kpi)  -- Chiara
 627     const Int_t ndght=2;
 628     nmasses=2;
 629     masses=new Float_t[nmasses];
 630     Int_t pdgdaughtersD0[ndght]={211,321};//pi,K
 631     masses[0]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0); //D0
 632     Int_t pdgdaughtersD0bar[ndght]={321,211};//K,pi
 633     masses[1]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0bar); //D0bar
 634   }
 635   if(fDecChannel==2){
 636     //D* -- Robert,Yifei, Alessandro
 637     nmasses=1;
 638     masses=new Float_t[nmasses];
 639     masses[0]=((AliAODRecoCascadeHF*)d)->DeltaInvMass();
 640   }
 641 }
 642
 643 //******************************************************************************
 644
 645 //Methods to fill the histograms, one for each channel
 646 //NB: the implementation for each candidate is responsibility of the corresponding developer
 647
 648 //******************************************************************************
 649 void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr, Double_t phiD){
 650   //D+ channel
 651   if(!isSel){
 652     if(fDebug>3)AliWarning("Candidate not selected\n");
 653     return;
 654   }
 655   if(!masses){
 656     if(fDebug>3)AliWarning("Masses not calculated\n");
 657     return;
 658   }
 659   Int_t icentrmin=icentr-fCentBinSizePerMil;
 660   ((TH2F*)fOutput->FindObject(Form("hMc2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 661   ((TH2F*)fOutput->FindObject(Form("hMs2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*deltaphi),masses[0]);
 662   ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
 663   ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentrmin,icentr)))->Fill(d->Pt(),masses[0]);
 664   ((TH2F*)fOutput->FindObject(Form("hCos2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*phiD),masses[0]);
 665   ((TH2F*)fOutput->FindObject(Form("hSin2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*phiD),masses[0]);
 666
 667
 668   Int_t pdgdaughters[3] = {211,321,211};
 669
 670   if(fReadMC){
 671     Int_t lab=-1;
 672     if(fUseAfterBurner){
 673       Bool_t isSignal=fAfterBurner->GetIsSignal();
 674       if(isSignal)lab=10;
 675     }else {
 676       lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
 677     }
 678     if(lab>=0){ //signal
 679       ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 680       ((TH2F*)fOutput->FindObject(Form("hMdeltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
 681     } else{ //background
 682       ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 683       ((TH2F*)fOutput->FindObject(Form("hMdeltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
 684     }
 685   }
 686 }
 687
 688 //******************************************************************************
 689 void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr,Double_t phiD){
 690
 691   //D0->Kpi channel
 692
 693   //mass histograms
 694   if(!masses){
 695     if(fDebug>3)AliWarning("Masses not calculated\n");
 696     return;
 697   }
 698   Int_t icentrmin=icentr-fCentBinSizePerMil;
 699   if(isSel==1 || isSel==3) {
 700     ((TH2F*)fOutput->FindObject(Form("hMc2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 701     ((TH2F*)fOutput->FindObject(Form("hMs2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*deltaphi),masses[0]);
 702     ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
 703     ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentrmin,icentr)))->Fill(d->Pt(),masses[0]);
 704     ((TH2F*)fOutput->FindObject(Form("hCos2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*phiD),masses[0]);
 705     ((TH2F*)fOutput->FindObject(Form("hSin2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*phiD),masses[0]);
 706   }
 707   if(isSel>=2) {
 708     ((TH2F*)fOutput->FindObject(Form("hMc2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 709     ((TH2F*)fOutput->FindObject(Form("hMs2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*deltaphi),masses[1]);
 710     ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[1]);
 711     ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentrmin,icentr)))->Fill(d->Pt(),masses[1]);
 712     ((TH2F*)fOutput->FindObject(Form("hCos2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*phiD),masses[1]);
 713     ((TH2F*)fOutput->FindObject(Form("hSin2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*phiD),masses[1]);
 714   }
 715
 716   //MC histograms
 717   if(fReadMC){
 718
 719     Int_t matchtoMC=-1;
 720
 721     //D0
 722     Int_t pdgdaughters[2];
 723     pdgdaughters[0]=211;//pi
 724     pdgdaughters[1]=321;//K
 725     Int_t nprongs=2;
 726     Int_t absPdgMom=421;
 727
 728     matchtoMC = d->MatchToMC(absPdgMom,arrayMC,nprongs,pdgdaughters);
 729
 730     Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
 731     if((isSel==1 || isSel==3)){ //D0
 732       if(matchtoMC>=0){
 733         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 734         Int_t pdgMC = dMC->GetPdgCode();
 735
 736         if(pdgMC==prongPdgPlus) {
 737           ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 738           ((TH2F*)fOutput->FindObject(Form("hMdeltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
 739         }
 740         else {
 741           ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiR_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 742           ((TH2F*)fOutput->FindObject(Form("hMdeltaphiR_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
 743         }
 744       } else {
 745         ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 746         ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
 747       }
 748     }
 749     if(isSel>=2){ //D0bar
 750       if(matchtoMC>=0){
 751         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 752         Int_t pdgMC = dMC->GetPdgCode();
 753
 754         if(pdgMC==prongPdgMinus) {
 755           ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 756           ((TH2F*)fOutput->FindObject(Form("hMdeltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[1]);
 757         }
 758         else {
 759           ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiR_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 760           ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[1]);
 761         }
 762       } else {
 763         ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 764         ((TH2F*)fOutput->FindObject(Form("hMdeltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[1]);
 765       }
 766     }
 767   }
 768 }
 769 //******************************************************************************
 770 void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr, Double_t phiD){
 771   //D* channel
 772   if(!isSel){
 773     if(fDebug>3)AliWarning("Candidate not selected\n");
 774     return;
 775   }
 776   if(!masses){
 777     if(fDebug>3)AliWarning("Masses not calculated\n");
 778     return;
 779   }
 780
 781   Float_t deltaphi1 = deltaphi;
 782   if(deltaphi1 > TMath::PiOver2()) deltaphi1 = TMath::Pi() - deltaphi1;
 783   Int_t icentrmin=icentr-fCentBinSizePerMil;
 784   ((TH2F*)fOutput->FindObject(Form("hMc2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 785   ((TH2F*)fOutput->FindObject(Form("hMs2deltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*deltaphi),masses[0]);
 786   ((TH2F*)fOutput->FindObject(Form("hMdeltaphi_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi1,masses[0]);
 787   ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentrmin,icentr)))->Fill(d->Pt(),masses[0]);
 788   ((TH2F*)fOutput->FindObject(Form("hCos2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*phiD),masses[0]);
 789   ((TH2F*)fOutput->FindObject(Form("hSin2phiMass_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Sin(2*phiD),masses[0]);
 790   Int_t pdgDgDStartoD0pi[2]={421,211};
 791   Int_t pdgDgD0toKpi[2]={321,211};
 792
 793   if(fReadMC){
 794     Int_t lab=-1;
 795     lab = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,arrayMC);
 796     if(lab>=0){ //signal
 797       ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 798       ((TH2F*)fOutput->FindObject(Form("hMdeltaphiS_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
 799     } else{ //background
 800       ((TH2F*)fOutput->FindObject(Form("hMc2deltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 801       ((TH2F*)fOutput->FindObject(Form("hMdeltaphiB_pt%dcentr%d_%d",ptbin,icentrmin,icentr)))->Fill(deltaphi,masses[0]);
 802     }
 803   }
 804 }
 805
 806 //________________________________________________________________________
 807 void AliAnalysisTaskSEHFv2::SetEventPlaneMethod(Int_t method){
 808   if(method>kVZEROC||method<0){
 809     AliWarning("No EP method associated to the selection, setting to TPC EP\n");
 810     method=kTPCVZERO;
 811   }
 812   fEventPlaneMeth=method;
 813 }
 814 //________________________________________________________________________
 815 void AliAnalysisTaskSEHFv2::SetNTPCSubEvents(Int_t nsub){
 816   if(nsub>3||nsub<2){
 817     AliWarning("Only 2 or 3 subevents implemented. Setting 2 subevents\n");
 818     nsub=2;
 819   }
 820   if(fEventPlaneMeth==kTPC&&nsub==3){
 821     AliWarning("V0 must be enabled to run 3 TPC subevents. Enabling...\n");
 822     fEventPlaneMeth=kTPCVZERO;
 823   }
 824   fSubEvents=nsub;
 825 }
 826 //________________________________________________________________________
 827 Float_t AliAnalysisTaskSEHFv2::GetPhi0Pi(Float_t phi){
 828   // Sets the phi angle in the range 0-pi
 829   Float_t result=phi;
 830   while(result<0){
 831     result=result+TMath::Pi();
 832   }
 833   while(result>TMath::Pi()){
 834     result=result-TMath::Pi();
 835   }
 836    return result;
 837 }
 838
 839 //________________________________________________________________________
 840 Float_t AliAnalysisTaskSEHFv2::GetEventPlaneForCandidate(AliAODRecoDecayHF* d, const TVector2* q,AliEventplane *pl, const TVector2* qsub1, const TVector2* qsub2){
 841   // remove autocorrelations
 842
 843   TArrayF* qx = 0x0;
 844   TArrayF* qy = 0x0;
 845   TVector2 qcopy;
 846   if(!fEtaGap){
 847     qx = pl->GetQContributionXArray();
 848     qy = pl->GetQContributionYArray();
 849     qcopy = *q;
 850   }else {
 851     if(d->Eta()<0.){
 852       qx = pl->GetQContributionXArraysub1();
 853       qy = pl->GetQContributionYArraysub1();
 854       qcopy = *qsub1;
 855     }else{
 856       qx = pl->GetQContributionXArraysub2();
 857       qy = pl->GetQContributionYArraysub2();
 858       qcopy = *qsub2;
 859     }
 860   }
 861
 862   if(fDecChannel==2){
 863     //D* -- Yifei, Alessandro,Robert
 864     AliAODRecoDecayHF2Prong* theD0particle = ((AliAODRecoCascadeHF*)d)->Get2Prong();
 865     AliAODTrack *track0 = (AliAODTrack*)theD0particle->GetDaughter(0);
 866     AliAODTrack *track1 = (AliAODTrack*)theD0particle->GetDaughter(1);
 867     AliAODTrack *track2 = ((AliAODRecoCascadeHF*)d)->GetBachelor();
 868     // reduce global q vector
 869
 870     TVector2 q0;
 871     if((track0->GetID()) < qx->fN){
 872       q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
 873
 874     TVector2 q1;
 875     if((track1->GetID()) < qx->fN){
 876       q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
 877
 878     TVector2 q2;
 879     if((track2->GetID()) < qx->fN){
 880       q2.Set(qx->At(track2->GetID()),qy->At(track2->GetID()));}
 881
 882     qcopy = qcopy -(q0+q1+q2);
 883
 884   }
 885
 886   // reduce Q vector for D+ and D0
 887
 888   if(fDecChannel==1){
 889     //D0 -- Chiara
 890     AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
 891     AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
 892
 893     TVector2 q0;
 894     if((track0->GetID()) < qx->fN){
 895       q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
 896
 897     TVector2 q1;
 898     if((track1->GetID()) < qx->fN){
 899       q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
 900
 901     qcopy = qcopy -(q0+q1);
 902   }
 903
 904   if(fDecChannel==0){
 905     //D+ -- Giacomo
 906     AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
 907     AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
 908     AliAODTrack *track2 = (AliAODTrack*)d->GetDaughter(2);
 909
 910     TVector2 q0;
 911     if((track0->GetID()) < qx->fN){
 912       q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
 913
 914     TVector2 q1;
 915     if((track1->GetID()) < qx->fN){
 916       q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
 917
 918     TVector2 q2;
 919     if((track2->GetID()) < qx->fN){
 920       q2.Set(qx->At(track2->GetID()),qy->At(track2->GetID()));}
 921
 922     qcopy = qcopy -(q0+q1+q2);
 923
 924   }
 925
 926   return qcopy.Phi()/2.;
 927
 928 }
 929 // //________________________________________________________________________
 930 // Float_t AliAnalysisTaskSEHFv2::GetEventPlaneFromV0(AliAODEvent *aodEvent){
 931 //   // Compute event plane for VZERO - Obsolete, used for 2010 data
 932
 933 //   Int_t centr=fRDCuts->GetCentrality(aodEvent);
 934 //   centr=centr-centr%10;
 935 //   //temporary fix
 936 //   if(centr<20)centr=20;
 937 //   if(centr>70)centr=70;
 938 //   //end temporary fix
 939 //   Int_t binx=0;
 940 //   Int_t iParHist=(centr-20)/10;
 941
 942 //   TString name;name.Form("parhist%d_%d",centr,centr+10);
 943
 944 //   if(fDebug>15)printf("EPfromV0 centr %d, iparhist %d (%p-%p)\n",centr,iParHist,fParHist->FindObject(name.Data()),fParHist->At(iParHist));
 945
 946 //   Int_t runnumber=aodEvent->GetRunNumber();
 947 //   if(fParHist->At(iParHist)){
 948 //     for(Int_t i=1;i<=((TH2D*)fParHist->At(iParHist))->GetNbinsX()&&binx<=0;i++){
 949 //       Int_t run=atoi(((TH2D*)fParHist->At(iParHist))->GetXaxis()->GetBinLabel(i));
 950 //       if(run>=runnumber)binx=i;
 951 //     }
 952 //   }else{
 953 //     fhEventsInfo->Fill(7);
 954 //   }
 955
 956 //   AliFlowTrackCuts* cutsRP = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts();
 957 //   cutsRP->SetEvent(aodEvent, MCEvent());//, 0x0);
 958 //   cutsRP->SetName( Form("rp_cuts") );
 959 //   AliFlowTrackCuts* dummy = new AliFlowTrackCuts("null_cuts");
 960 //   dummy->SetParamType(AliFlowTrackCuts::kGlobal);
 961 //   dummy->SetPtRange(+1,-1); // select nothing QUICK
 962 //   dummy->SetEtaRange(+1,-1); // select nothing VZERO
 963 //   dummy->SetEvent(aodEvent,MCEvent());
 964
 965 //   //////////////// construct the flow event container ////////////
 966 //   AliFlowEvent flowEvent(cutsRP,dummy);
 967 //   flowEvent.SetReferenceMultiplicity( 64 );
 968 //   for(Int_t i=0;i<64&&binx>0;i++){
 969 //     AliFlowTrack *flowTrack=flowEvent.GetTrack(i);
 970 //     Double_t inte=((TH2D*)fParHist->At(iParHist))->Integral(binx,binx,i+1,i+1);
 971 //     if(inte>0)flowTrack->SetWeight(flowTrack->Weight()/inte);
 972 //   }
 973 //   if(fDebug>15)printf("EPfromV0 flow tracks weights done\n");
 974
 975 //   AliFlowVector qvec=flowEvent.GetQ(fV0EPorder);
 976 //   Double_t angleEP=(1./(Double_t)fV0EPorder)*qvec.Phi();
 977 //   if(fDebug>15)printf("EPfromV0 phi %f\n",angleEP);
 978 //   return angleEP;
 979 // }
 980 //________________________________________________________________________
 981 void AliAnalysisTaskSEHFv2::Terminate(Option_t */*option*/)
 982 {
 983   // Terminate analysis
 984   //
 985   if(fDebug > 1) printf("AnalysisTaskSEHFv2: Terminate() \n");
 986
 987   return;
 988 }
 989