PWG3/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   fParHist(0),
  87   fLowmasslimit(1.765),
  88   fUpmasslimit(1.965),
  89   fNPtBins(1),
  90   fNPhiBinLims(2),
  91   fPhiBins(0),
  92   fNMassBins(200),
  93   fReadMC(kFALSE),
  94   fUseAfterBurner(kFALSE),
  95   fDecChannel(0),
  96   fAfterBurner(0),
  97   fUseV0EP(kFALSE),
  98   fV0EPorder(2),
  99   fMinCentr(20),
 100   fMaxCentr(80),
 101   fEtaGap(kFALSE)
 102 {
 103   // Default constructor
 104 }
 105
 106 //________________________________________________________________________
 107 AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2(const char *name,AliRDHFCuts *rdCuts,Int_t decaychannel, Int_t nlimsphibin, Float_t *phibinlimits):
 108   AliAnalysisTaskSE(name),
 109   fhEventsInfo(0),
 110   fOutput(0),
 111   fRDCuts(rdCuts),
 112   fParHist(0),
 113   fLowmasslimit(0),
 114   fUpmasslimit(0),
 115   fNPtBins(1),
 116   fNPhiBinLims(2),
 117   fPhiBins(0),
 118   fNMassBins(200),
 119   fReadMC(kFALSE),
 120   fUseAfterBurner(kFALSE),
 121   fDecChannel(decaychannel),
 122   fAfterBurner(0),
 123   fUseV0EP(kFALSE),
 124   fV0EPorder(2),
 125   fMinCentr(20),
 126   fMaxCentr(80),
 127   fEtaGap(kFALSE)
 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.135,(Float_t)0.165);
 153   else SetMassLimits((Float_t)0.2,pdg); //check range
 154   fNPtBins=fRDCuts->GetNPtBins();
 155   if(nlimsphibin>2) fNPhiBinLims=nlimsphibin;
 156   else AliInfo("At least 2 limits in Delta phi needed");
 157
 158   fPhiBins=new Float_t[fNPhiBinLims];
 159   for(Int_t i=0;i<fNPhiBinLims;i++) fPhiBins[i]=phibinlimits[i];
 160
 161   if(fDebug>1)fRDCuts->PrintAll();
 162   // Output slot #1 writes into a TH1F container
 163   DefineOutput(1,TH1F::Class());   //Info on the number of events etc.
 164   // Output slot #2 writes into a TList container
 165   DefineOutput(2,TList::Class());  //Main output
 166   // Output slot #3 writes into a AliRDHFCuts container (cuts)
 167   switch(fDecChannel){
 168   case 0:
 169     DefineOutput(3,AliRDHFCutsDplustoKpipi::Class());  //Cut object for Dplus
 170     break;
 171   case 1:
 172     DefineOutput(3,AliRDHFCutsD0toKpi::Class());  //Cut object for D0
 173     break;
 174   case 2:
 175     DefineOutput(3,AliRDHFCutsDStartoKpipi::Class());  //Cut object for D*
 176     break;
 177   }
 178   //DefineOutput(4,AliFlowEventSimple::Class());
 179   //DefineOutput(4,TList::Class());
 180 }
 181
 182 //________________________________________________________________________
 183 AliAnalysisTaskSEHFv2::~AliAnalysisTaskSEHFv2()
 184 {
 185   // Destructor
 186   delete fOutput;
 187   delete fhEventsInfo;
 188   delete fRDCuts;
 189   delete fParHist;
 190
 191   for(Int_t i=0;i<6;i++){
 192     delete fHistvzero[i];
 193   }
 194
 195   delete fAfterBurner;
 196 }
 197 //_________________________________________________________________
 198 void  AliAnalysisTaskSEHFv2::SetVZEROParHist(TH2D** histPar){
 199   // Set the histograms for VZERO EP parameters
 200   for(Int_t i=0;i<6;i++)fHistvzero[i]=(TH2D*)histPar[i]->Clone();
 201   for(Int_t i=0;i<6;i++){
 202     if(!fHistvzero[i]){
 203       printf("No VZERO histograms!\n");
 204       fUseV0EP=kFALSE;
 205       return;
 206     }
 207   }
 208   DefineOutput(4,TList::Class());
 209   fUseV0EP=kTRUE;
 210 }
 211 //_________________________________________________________________
 212 void  AliAnalysisTaskSEHFv2::SetMassLimits(Float_t range, Int_t pdg){
 213   // Set limits for mass spectra plots
 214   Float_t mass=0;
 215   Int_t abspdg=TMath::Abs(pdg);
 216   mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
 217   fUpmasslimit = mass+range;
 218   fLowmasslimit = mass-range;
 219 }
 220 //_________________________________________________________________
 221 void  AliAnalysisTaskSEHFv2::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 222   // Set limits for mass spectra plots
 223   if(uplimit>lowlimit)
 224     {
 225       fUpmasslimit = uplimit;
 226       fLowmasslimit = lowlimit;
 227     }
 228 }
 229
 230
 231 //________________________________________________________________________
 232 void AliAnalysisTaskSEHFv2::LocalInit()
 233 {
 234   // Initialization
 235
 236   if(fDebug > 1) printf("AnalysisTaskSEHFv2::Init() \n");
 237
 238   fRDCuts->SetMinCentrality(fMinCentr);
 239   fRDCuts->SetMaxCentrality(fMaxCentr);
 240
 241   switch(fDecChannel){
 242   case 0:
 243     {
 244       AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCuts)));
 245       // Post the data
 246       PostData(3,copycut);
 247     }
 248     break;
 249   case 1:
 250     {
 251       AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCuts)));
 252       // Post the data
 253       PostData(3,copycut);
 254     }
 255     break;
 256   case 2:
 257     {
 258       AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCuts)));
 259       // Post the data
 260       PostData(3,copycut);
 261     }
 262     break;
 263   default:
 264     return;
 265   }
 266   return;
 267 }
 268 //________________________________________________________________________
 269 void AliAnalysisTaskSEHFv2::UserCreateOutputObjects()
 270 {
 271   // Create the output container
 272
 273   if(fDebug > 1) printf("AnalysisTaskSEHFv2::UserCreateOutputObjects() \n");
 274
 275   fhEventsInfo=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Number of AODs scanned",8,-0.5,7.5);
 276   fhEventsInfo->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 277   fhEventsInfo->GetXaxis()->SetBinLabel(2,"nEvSelected");
 278   fhEventsInfo->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
 279   fhEventsInfo->GetXaxis()->SetBinLabel(4,"out of pt bounds");
 280   fhEventsInfo->GetXaxis()->SetBinLabel(5,"Pile-up Rej");
 281   fhEventsInfo->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
 282   fhEventsInfo->GetXaxis()->SetBinLabel(7,Form("Ev Sel in Centr %.0f-%.0f%s",fRDCuts->GetMinCentrality(),fRDCuts->GetMaxCentrality(),"%"));
 283   fhEventsInfo->GetXaxis()->SetBinLabel(8,"mismatch lab");
 284   fhEventsInfo->GetXaxis()->SetNdivisions(1,kFALSE);
 285
 286
 287   // Several histograms are more conveniently managed in a TList
 288   fOutput = new TList();
 289   fOutput->SetOwner();
 290   fOutput->SetName("MainOutput");
 291
 292   for(Int_t icentr=fMinCentr+5;icentr<=fMaxCentr;icentr=icentr+5){
 293     TString centrname;centrname.Form("centr%d_%d",icentr-5,icentr);
 294     for(Int_t i=0;i<fNPtBins;i++){
 295
 296       TString hname;
 297       TString title;
 298       hname.Form("hPhi_pt%d",i);hname.Append(centrname.Data());
 299       title.Form("Phi distribution (Pt bin %d %s);#phi;Entries",i,centrname.Data());
 300
 301       TH1F* hPhi=new TH1F(hname.Data(),title.Data(),96,0.,2*TMath::Pi());
 302       hPhi->Sumw2();
 303       fOutput->Add(hPhi);
 304
 305       for(Int_t j=0;j<fNPhiBinLims-1;j++){
 306
 307         hname.Form("hMass_pt%dphi%d",i,j);hname.Append(centrname.Data());
 308         title.Form("Invariant mass (Pt bin %d, Phi bin %d, %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
 309
 310         TH1F* hMass=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
 311         hMass->Sumw2();
 312         fOutput->Add(hMass);
 313
 314
 315         if(fReadMC){
 316           hname.Form("hSgn_pt%dphi%d",i,j);hname.Append(centrname.Data());
 317           title.Form("Invariant mass S (Pt bin %d, Phi bin %d, %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
 318           TH1F* hSgn=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
 319           hSgn->Sumw2();
 320           fOutput->Add(hSgn);
 321
 322           hname.Form("hBkg_pt%dphi%d",i,j);hname.Append(centrname.Data());
 323           title.Form("Invariant mass B (Pt bin %d, Phi bin %d, %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
 324           TH1F* hBkg=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
 325           hBkg->Sumw2();
 326           fOutput->Add(hBkg);
 327
 328           if((fDecChannel != AliAnalysisTaskSEHFv2::kDplustoKpipi) && (fDecChannel != AliAnalysisTaskSEHFv2::kDstartoKpipi)){
 329             hname.Form("hRfl_pt%dphi%d",i,j);hname.Append(centrname.Data());
 330             title.Form("Invariant mass Reflections (Pt bin %d, Phi bin %d %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
 331             TH1F* hRfl=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
 332             hRfl->Sumw2();
 333             fOutput->Add(hRfl);
 334           }
 335         }
 336       }
 337
 338       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);
 339       fOutput->Add(hMc2phi);
 340
 341       if (fReadMC){
 342         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);
 343         fOutput->Add(hMc2phiS);
 344         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);
 345         fOutput->Add(hMphiS);
 346         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);
 347         fOutput->Add(hMc2phiB);
 348         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);
 349         fOutput->Add(hMphiB);
 350         if((fDecChannel != AliAnalysisTaskSEHFv2::kDplustoKpipi) &&(fDecChannel != AliAnalysisTaskSEHFv2::kDstartoKpipi)){
 351           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);
 352           fOutput->Add(hMc2phiR);
 353           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);
 354           fOutput->Add(hMphiR);
 355         }
 356       }
 357     }
 358
 359     TH2F* hMphi=new TH2F(Form("hMphi%s",centrname.Data()),Form("Mass vs #Delta#phi %s;#Delta#phi;M (GeV/c^{2})",centrname.Data()),96,0,TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
 360     fOutput->Add(hMphi);
 361
 362     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);
 363     fOutput->Add(hMPtCand);
 364
 365     TH1F* hEvPlaneneg=new TH1F(Form("hEvPlaneneg%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
 366     fOutput->Add(hEvPlaneneg);
 367
 368         TH1F* hEvPlanepos=new TH1F(Form("hEvPlanepos%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
 369     fOutput->Add(hEvPlanepos);
 370
 371     //TH1F* hEvPlaneCheck=new TH1F(Form("hEvPlaneCheck%s",centrname.Data()),Form("Event plane angle - Event plane angle per candidate %s;(#phi(Ev Plane) - #phi(Ev Plane Candidate))/#phi(EvPlane);Entries",centrname.Data()),200,-0.2,0.2);
 372     //fOutput->Add(hEvPlaneCheck);
 373
 374     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());
 375     fOutput->Add(hEvPlaneCand);
 376
 377     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);
 378     fOutput->Add(hEvPlaneReso);
 379   }
 380
 381   TH1F* hPhiBins=new TH1F("hPhiBins","Bins in #Delta#phi used in this analysis;#phi bin;n jobs",fNPhiBinLims-1,fPhiBins);
 382   fOutput->Add(hPhiBins);
 383   for(Int_t k=0;k<fNPhiBinLims-1;k++)hPhiBins->SetBinContent(k+1,1);
 384   PostData(1,fhEventsInfo);
 385   PostData(2,fOutput);
 386   if(fUseV0EP){
 387     fParHist = new TList();
 388     fParHist->SetOwner();
 389     fParHist->SetName("VZEROcorr");
 390     for(Int_t i=0;i<6;i++){
 391       fParHist->Add((TH2D*)fHistvzero[i]);
 392     }
 393     PostData(4,fParHist);
 394   }
 395   return;
 396 }
 397
 398 //________________________________________________________________________
 399 void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
 400 {
 401   // Execute analysis for current event:
 402   // heavy flavor candidates association to MC truth
 403   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 404   if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
 405   // Post the data already here
 406   PostData(1,fhEventsInfo);
 407   PostData(2,fOutput);
 408
 409   TClonesArray *arrayProng =0;
 410   Int_t absPdgMom=0;
 411   if(!aod && AODEvent() && IsStandardAOD()) {
 412     // In case there is an AOD handler writing a standard AOD, use the AOD
 413     // event in memory rather than the input (ESD) event.
 414     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 415     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 416     // have to taken from the AOD event hold by the AliAODExtension
 417     AliAODHandler* aodHandler = (AliAODHandler*)
 418       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 419     if(aodHandler->GetExtensions()) {
 420
 421       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 422       AliAODEvent *aodFromExt = ext->GetAOD();
 423
 424
 425       if(fDecChannel==0){
 426         absPdgMom=411;
 427         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 428       }
 429       if(fDecChannel==1){
 430         absPdgMom=421;
 431         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 432       }
 433       if(fDecChannel==2){
 434         absPdgMom=413;
 435         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 436       }
 437     }
 438   } else {
 439     if(fDecChannel==0){
 440       absPdgMom=411;
 441       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 442     }
 443     if(fDecChannel==1){
 444       absPdgMom=421;
 445       arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 446     }
 447     if(fDecChannel==2){
 448       absPdgMom=413;
 449       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
 450     }
 451   }
 452
 453   if(!arrayProng) {
 454     AliError("AliAnalysisTaskSEHFv2::UserExec:Branch not found!\n");
 455     return;
 456   }
 457
 458   // fix for temporary bug in ESDfilter
 459   // the AODs with null vertex pointer didn't pass the PhysSel
 460   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 461
 462   TClonesArray *arrayMC=0;
 463   AliAODMCHeader *mcHeader=0;
 464
 465   // load MC particles
 466   if(fReadMC){
 467
 468     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 469     if(!arrayMC) {
 470       AliWarning("AliAnalysisTaskSEHFv2::UserExec:MC particles branch not found!\n");
 471       //    return;
 472     }
 473
 474     // load MC header
 475     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 476     if(!mcHeader) {
 477       AliError("AliAnalysisTaskSEHFv2::UserExec:MC header branch not found!\n");
 478       return;
 479     }
 480   }
 481
 482   fhEventsInfo->Fill(0); // count event
 483
 484   AliAODRecoDecayHF *d=0;
 485
 486   Int_t nCand = arrayProng->GetEntriesFast();
 487   if(fDebug>2) printf("Number of D2H: %d\n",nCand);
 488
 489   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
 490   TString trigclass=aod->GetFiredTriggerClasses();
 491   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fhEventsInfo->Fill(5);
 492
 493   if(fRDCuts->IsEventSelectedInCentrality(aod)>0) return;
 494   else fhEventsInfo->Fill(6);
 495   if(fRDCuts->IsEventSelected(aod)) fhEventsInfo->Fill(1);
 496   else{
 497     if(fRDCuts->GetWhyRejection()==1) // rejected for pileup
 498       fhEventsInfo->Fill(4);
 499     return;
 500   }
 501
 502   AliEventplane *pl=0x0;
 503   TVector2* q=0x0;
 504   Double_t rpangleevent=0;
 505   Double_t rpangleeventneg=0;
 506   Double_t rpangleeventpos=0;
 507   Double_t eventplane=0;
 508   TVector2 *qsub1=0x0;
 509   TVector2 *qsub2=0x0;
 510
 511   //determine centrality bin
 512   Float_t centr=fRDCuts->GetCentrality(aod);
 513   Int_t icentr=0;
 514   for(Int_t ic=fMinCentr+5;ic<=fMaxCentr;ic=ic+5){
 515     if(ic>centr){
 516       icentr=ic;
 517       break;
 518     }
 519   }
 520   TString centrbinname=Form("centr%d_%d",icentr-5,icentr);
 521   if(fReadMC){
 522     fUseV0EP=kTRUE;
 523     TRandom3 *g = new TRandom3(0);
 524     rpangleevent=g->Rndm()*TMath::Pi();
 525     delete g;g=0x0;
 526     eventplane=rpangleevent;
 527         ((TH1F*)fOutput->FindObject(Form("hEvPlanepos%s",centrbinname.Data())))->Fill(rpangleevent);
 528     if(fUseAfterBurner)fAfterBurner->SetEventPlane((Double_t)rpangleevent);
 529   }else{
 530     if(fUseV0EP){
 531       rpangleevent=GetEventPlaneFromV0(aod);
 532       eventplane=rpangleevent;
 533           ((TH1F*)fOutput->FindObject(Form("hEvPlanepos%s",centrbinname.Data())))->Fill(rpangleevent);
 534     }else{
 535       // event plane and resolution
 536       //--------------------------------------------------------------------------
 537       // extracting Q vectors and calculating v2 + resolution
 538       pl = aod->GetHeader()->GetEventplaneP();
 539       if(!pl){
 540         AliError("AliAnalysisTaskSEHFv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
 541         return;
 542       }
 543       if(fEtaGap){
 544         qsub1 = pl->GetQsub1();
 545         qsub2 = pl->GetQsub2();
 546         rpangleeventpos = qsub1->Phi()/2.;
 547         rpangleeventneg = qsub2->Phi()/2.;
 548         ((TH1F*)fOutput->FindObject(Form("hEvPlanepos%s",centrbinname.Data())))->Fill(rpangleeventpos);
 549         ((TH1F*)fOutput->FindObject(Form("hEvPlaneneg%s",centrbinname.Data())))->Fill(rpangleeventneg);
 550         }
 551       else if(!fEtaGap){
 552         q = pl->GetQVector();
 553         rpangleevent = pl->GetEventplane("Q");
 554         ((TH1F*)fOutput->FindObject(Form("hEvPlanepos%s",centrbinname.Data())))->Fill(rpangleevent); // reaction plane angle without autocorrelations removal
 555         }
 556
 557       Double_t deltaPsi = pl->GetQsubRes();
 558       if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){
 559         if(deltaPsi>0.) deltaPsi-=TMath::Pi();
 560         else deltaPsi +=TMath::Pi();
 561       } // difference of subevents reaction plane angle cannot be bigger than phi/2
 562       Double_t planereso = TMath::Cos(2.*deltaPsi); // reaction plane resolution
 563       //--------------------------------------------------------------------------
 564       ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso%s",centrbinname.Data())))->Fill(planereso);
 565     }
 566   }
 567
 568
 569   for (Int_t iCand = 0; iCand < nCand; iCand++) {
 570
 571     d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iCand);
 572
 573     Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(absPdgMom));
 574     Int_t isSelected= fRDCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 575     Bool_t isSelBit=kTRUE;
 576     if(fDecChannel==0) isSelBit=d->HasSelectionBit(AliRDHFCuts::kDplusCuts);
 577     if(fDecChannel==1) isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts);
 578     if(fDecChannel==2) isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
 579     if(isSelected && isFidAcc && isSelBit) {
 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       Int_t ptbin=fRDCuts->PtBin(d->Pt());
 588       if(ptbin==-1) {
 589         fhEventsInfo->Fill(3);
 590         continue;
 591       }
 592
 593       if(!fUseV0EP) {
 594         eventplane = GetEventPlaneForCandidate(d,q,pl,qsub1,qsub2); // remove autocorrelations
 595         ((TH1F*)fOutput->FindObject(Form("hEvPlaneCand%s",centrbinname.Data())))->Fill(rpangleevent-eventplane);
 596         //((TH1F*)fOutput->FindObject(Form("hEvPlaneCheck%s",centrbinname.Data())))->Fill((rpangleevent-eventplane)/100.*rpangleevent);
 597       }
 598
 599       Float_t phi=d->Phi();
 600       ((TH1F*)fOutput->FindObject(Form("hPhi_pt%d%s",ptbin,centrbinname.Data())))->Fill(phi);
 601
 602       if(fReadMC&&fUseAfterBurner)phi=fAfterBurner->GetNewAngle(d,arrayMC);
 603       Float_t deltaphi=GetPhi0Pi(phi-eventplane);
 604
 605       //fill the histograms with the appropriate method
 606       if(fDecChannel==0)FillDplus(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
 607       if(fDecChannel==1)FillD02p(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
 608       if(fDecChannel==2)FillDstar(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
 609     }// end if selected
 610   }
 611
 612   PostData(1,fhEventsInfo);
 613   PostData(2,fOutput);
 614   //PostData(4,flowEvent);
 615   return;
 616 }
 617
 618 //***************************************************************************
 619
 620 // Methods used in the UserExec
 621
 622 void AliAnalysisTaskSEHFv2::CalculateInvMasses(AliAODRecoDecayHF* d,Float_t*& masses,Int_t& nmasses){
 623   //Calculates all the possible invariant masses for each candidate
 624   //NB: the implementation for each candidate is responsibility of the corresponding developer
 625
 626   if(fDecChannel==0){
 627     //D+ -- Giacomo
 628     nmasses=1;
 629     masses=new Float_t[nmasses];
 630     Int_t pdgdaughters[3] = {211,321,211};
 631     masses[0]=d->InvMass(3,(UInt_t*)pdgdaughters);
 632   }
 633   if(fDecChannel==1){
 634     //D0 (Kpi)  -- Chiara
 635     const Int_t ndght=2;
 636     nmasses=2;
 637     masses=new Float_t[nmasses];
 638     Int_t pdgdaughtersD0[ndght]={211,321};//pi,K
 639     masses[0]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0); //D0
 640     Int_t pdgdaughtersD0bar[ndght]={321,211};//K,pi
 641     masses[1]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0bar); //D0bar
 642   }
 643   if(fDecChannel==2){
 644     //D* -- Robert,Yifei, Alessandro
 645     nmasses=1;
 646     masses=new Float_t[nmasses];
 647     masses[0]=((AliAODRecoCascadeHF*)d)->DeltaInvMass();
 648   }
 649 }
 650
 651 //******************************************************************************
 652
 653 //Methods to fill the istograms with MC information, one for each candidate
 654 //NB: the implementation for each candidate is responsibility of the corresponding developer
 655
 656 void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
 657   //D+ channel
 658   if(!isSel){
 659     if(fDebug>3)AliWarning("Candidate not selected\n");
 660     return;
 661   }
 662   if(!masses){
 663     if(fDebug>3)AliWarning("Masses not calculated\n");
 664     return;
 665   }
 666
 667   Int_t phibin=GetPhiBin(deltaphi);
 668
 669   ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
 670   ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 671   ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 672   ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
 673   Int_t pdgdaughters[3] = {211,321,211};
 674
 675   if(fReadMC){
 676     Int_t lab=-1;
 677     if(fUseAfterBurner){
 678       Bool_t isSignal=fAfterBurner->GetIsSignal();
 679       if(isSignal)lab=10;
 680     }else {
 681       lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
 682     }
 683     if(lab>=0){ //signal
 684       ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
 685       ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 686       ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 687     } else{ //background
 688       ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
 689       ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 690       ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 691     }
 692   }
 693 }
 694
 695
 696 void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
 697
 698   //D0->Kpi channel
 699
 700   //mass histograms
 701   if(!masses){
 702     if(fDebug>3)AliWarning("Masses not calculated\n");
 703     return;
 704   }
 705   Int_t phibin=GetPhiBin(deltaphi);
 706   if(isSel==1 || isSel==3) {
 707     ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
 708     ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 709     ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 710     ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
 711   }
 712   if(isSel>=2) {
 713     ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);
 714     ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 715     ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[1]);
 716     ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[1]);
 717   }
 718
 719
 720   //MC histograms
 721   if(fReadMC){
 722
 723     Int_t matchtoMC=-1;
 724
 725     //D0
 726     Int_t pdgdaughters[2];
 727     pdgdaughters[0]=211;//pi
 728     pdgdaughters[1]=321;//K
 729     Int_t nprongs=2;
 730     Int_t absPdgMom=421;
 731
 732     matchtoMC = d->MatchToMC(absPdgMom,arrayMC,nprongs,pdgdaughters);
 733
 734     Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
 735     if((isSel==1 || isSel==3)){ //D0
 736       if(matchtoMC>=0){
 737         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 738         Int_t pdgMC = dMC->GetPdgCode();
 739
 740         if(pdgMC==prongPdgPlus) {
 741           ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
 742           ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 743           ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 744         }
 745         else {
 746           ((TH1F*)fOutput->FindObject(Form("hRfl_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
 747           ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 748           ((TH2F*)fOutput->FindObject(Form("hMphiRcentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 749         }
 750       } else {
 751         ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
 752         ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 753         ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 754       }
 755     }
 756     if(isSel>=2){ //D0bar
 757       if(matchtoMC>=0){
 758         AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
 759         Int_t pdgMC = dMC->GetPdgCode();
 760
 761         if(pdgMC==prongPdgMinus) {
 762           ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);
 763           ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 764           ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
 765         }
 766         else {
 767           ((TH1F*)fOutput->FindObject(Form("hRfl_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);
 768           ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 769           ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[1]);
 770         }
 771       } else {
 772         ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);
 773         ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
 774         ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
 775       }
 776     }
 777   }
 778 }
 779
 780 void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, const Float_t* masses,Int_t isSel,Int_t icentr){
 781   //D* channel
 782   if(!isSel){
 783     if(fDebug>3)AliWarning("Candidate not selected\n");
 784     return;
 785   }
 786   if(!masses){
 787     if(fDebug>3)AliWarning("Masses not calculated\n");
 788     return;
 789   }
 790   Int_t phibin=GetPhiBin(deltaphi);
 791
 792   ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
 793   ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 794   ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 795   ((TH2F*)fOutput->FindObject(Form("hMPtCandcentr%d_%d",icentr-5,icentr)))->Fill(d->Pt(),masses[0]);
 796   Int_t pdgDgDStartoD0pi[2]={421,211};
 797   Int_t pdgDgD0toKpi[2]={321,211};
 798
 799   if(fReadMC){
 800     Int_t lab=-1;
 801     lab = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,arrayMC);
 802     if(lab>=0){ //signal
 803       ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
 804       ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 805       ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 806     } else{ //background
 807       ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
 808       ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
 809       ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
 810     }
 811   }
 812
 813
 814 }
 815
 816
 817 //________________________________________________________________________
 818 Int_t AliAnalysisTaskSEHFv2::GetPhiBin(Float_t deltaphi) const{
 819
 820   //give the bin corresponding to the value of deltaphi according to the binning requested in the constructor
 821   Int_t phibin=0;
 822   for(Int_t i=0;i<fNPhiBinLims-1;i++) {
 823     if(deltaphi>=fPhiBins[i] && deltaphi<fPhiBins[i+1]) {
 824       phibin=i;
 825       break;
 826     }
 827   }
 828   return phibin;
 829 }
 830 //________________________________________________________________________
 831 // Float_t AliAnalysisTaskSEHFv2::GetPhi02Pi(Float_t phi){
 832 //   Float_t result=phi;
 833 //   while(result<0){
 834 //     result=result+2*TMath::Pi();
 835 //   }
 836 //   while(result>TMath::Pi()*2){
 837 //     result=result-2*TMath::Pi();
 838 //   }
 839 //   return result;
 840 // }
 841 //________________________________________________________________________
 842 Float_t AliAnalysisTaskSEHFv2::GetPhi0Pi(Float_t phi){
 843   // Sets the phi angle in the range 0-pi
 844   Float_t result=phi;
 845   while(result<0){
 846     result=result+TMath::Pi();
 847   }
 848   while(result>TMath::Pi()){
 849     result=result-TMath::Pi();
 850   }
 851   return result;
 852 }
 853
 854 //________________________________________________________________________
 855 Float_t AliAnalysisTaskSEHFv2::GetEventPlaneForCandidate(AliAODRecoDecayHF* d, const TVector2* q,AliEventplane *pl, const TVector2* qsub1, const TVector2* qsub2){
 856   // remove autocorrelations
 857
 858   TArrayF* qx = 0x0;
 859   TArrayF* qy = 0x0;
 860   TVector2 qcopy;
 861   if(!fEtaGap){
 862     qx = pl->GetQContributionXArray();
 863     qy = pl->GetQContributionYArray();
 864     qcopy = *q;
 865     }
 866   else {
 867     if(d->Eta()>0.){
 868       qx = pl->GetQContributionXArraysub1();
 869       qy = pl->GetQContributionYArraysub1();
 870       qcopy = *qsub1;
 871     }
 872     else{
 873       qx = pl->GetQContributionXArraysub2();
 874       qy = pl->GetQContributionYArraysub2();
 875       qcopy = *qsub2;
 876     }
 877   }
 878
 879
 880
 881   if(fDecChannel==2){
 882     //D* -- Yifei, Alessandro,Robert
 883     AliAODRecoDecayHF2Prong* theD0particle = ((AliAODRecoCascadeHF*)d)->Get2Prong();
 884     AliAODTrack *track0 = (AliAODTrack*)theD0particle->GetDaughter(0);
 885     AliAODTrack *track1 = (AliAODTrack*)theD0particle->GetDaughter(1);
 886     AliAODTrack *track2 = ((AliAODRecoCascadeHF*)d)->GetBachelor();
 887     // reduce global q vector
 888
 889     TVector2 q0;
 890     if((track0->GetID()) < qx->fN){
 891       q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
 892
 893     TVector2 q1;
 894     if((track1->GetID()) < qx->fN){
 895       q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
 896
 897     TVector2 q2;
 898     if((track2->GetID()) < qx->fN){
 899       q2.Set(qx->At(track2->GetID()),qy->At(track2->GetID()));}
 900
 901     qcopy = qcopy -(q0+q1+q2);
 902
 903   }
 904
 905   // reduce Q vector for D+ and D0
 906
 907   if(fDecChannel==1){
 908     //D0 -- Chiara
 909     AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
 910     AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
 911
 912     TVector2 q0;
 913     if((track0->GetID()) < qx->fN){
 914       q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
 915
 916     TVector2 q1;
 917     if((track1->GetID()) < qx->fN){
 918       q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
 919
 920     qcopy = qcopy -(q0+q1);
 921   }
 922
 923   if(fDecChannel==0){
 924     //D+ -- Giacomo
 925     AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
 926     AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
 927     AliAODTrack *track2 = (AliAODTrack*)d->GetDaughter(2);
 928
 929     TVector2 q0;
 930     if((track0->GetID()) < qx->fN){
 931       q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
 932
 933     TVector2 q1;
 934     if((track1->GetID()) < qx->fN){
 935       q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
 936
 937     TVector2 q2;
 938     if((track2->GetID()) < qx->fN){
 939       q2.Set(qx->At(track2->GetID()),qy->At(track2->GetID()));}
 940
 941     qcopy = qcopy -(q0+q1+q2);
 942
 943   }
 944
 945   return qcopy.Phi()/2.;
 946
 947 }
 948 //________________________________________________________________________
 949 Float_t AliAnalysisTaskSEHFv2::GetEventPlaneFromV0(AliAODEvent *aodEvent){
 950   // Compute event plane for VZERO
 951
 952   Int_t centr=fRDCuts->GetCentrality(aodEvent);
 953   centr=centr-centr%10;
 954   //temporary fix
 955   if(centr<20)centr=20;
 956   if(centr>70)centr=70;
 957   //end temporary fix
 958   Int_t binx=0;
 959   Int_t iParHist=(centr-20)/10;
 960
 961   TString name;name.Form("parhist%d_%d",centr,centr+10);
 962
 963   if(fDebug>15)printf("EPfromV0 centr %d, iparhist %d (%p-%p)\n",centr,iParHist,fParHist->FindObject(name.Data()),fParHist->At(iParHist));
 964
 965   Int_t runnumber=aodEvent->GetRunNumber();
 966   if(fParHist->At(iParHist)){
 967     for(Int_t i=1;i<=((TH2D*)fParHist->At(iParHist))->GetNbinsX()&&binx<=0;i++){
 968       Int_t run=atoi(((TH2D*)fParHist->At(iParHist))->GetXaxis()->GetBinLabel(i));
 969       if(run>=runnumber)binx=i;
 970     }
 971   }else{
 972     fhEventsInfo->Fill(7);
 973   }
 974
 975   AliFlowTrackCuts* cutsRP = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts();
 976   cutsRP->SetEvent(aodEvent, MCEvent());//, 0x0);
 977   cutsRP->SetName( Form("rp_cuts") );
 978   AliFlowTrackCuts* dummy = new AliFlowTrackCuts("null_cuts");
 979   dummy->SetParamType(AliFlowTrackCuts::kGlobal);
 980   dummy->SetPtRange(+1,-1); // select nothing QUICK
 981   dummy->SetEtaRange(+1,-1); // select nothing VZERO
 982   dummy->SetEvent(aodEvent,MCEvent());
 983
 984   //////////////// construct the flow event container ////////////
 985   AliFlowEvent flowEvent(cutsRP,dummy);
 986   flowEvent.SetReferenceMultiplicity( 64 );
 987   for(Int_t i=0;i<64&&binx>0;i++){
 988     AliFlowTrack *flowTrack=flowEvent.GetTrack(i);
 989     Double_t inte=((TH2D*)fParHist->At(iParHist))->Integral(binx,binx,i+1,i+1);
 990     if(inte>0)flowTrack->SetWeight(flowTrack->Weight()/inte);
 991   }
 992   if(fDebug>15)printf("EPfromV0 flow tracks weights done\n");
 993
 994   AliFlowVector qvec=flowEvent.GetQ(fV0EPorder);
 995   Double_t angleEP=(1./(Double_t)fV0EPorder)*qvec.Phi();
 996   if(fDebug>15)printf("EPfromV0 phi %f\n",angleEP);
 997   return angleEP;
 998 }
 999 //________________________________________________________________________
1000 void AliAnalysisTaskSEHFv2::Terminate(Option_t */*option*/)
1001 {
1002
1003   // Terminate analysis
1004   //
1005   if(fDebug > 1) printf("AnalysisTaskSEHFv2: Terminate() \n");
1006   /*
1007     fhEventsInfo = dynamic_cast<TH1F*> (GetOutputData(1));
1008     if(!fhEventsInfo){
1009     printf("Error: hEventsInfo not available\n");
1010     return;
1011     }
1012
1013     fOutput = dynamic_cast<TList*> (GetOutputData(2));
1014     if (!fOutput) {
1015     printf("ERROR: fOutput not available\n");
1016     return;
1017     }
1018     switch(fDecChannel){
1019     case 0:
1020     fRDCuts = dynamic_cast<AliRDHFCutsDplustoKpipi*> (GetOutputData(3));
1021     break;
1022     case 1:
1023     fRDCuts = dynamic_cast<AliRDHFCutsD0toKpi*> (GetOutputData(3));
1024     break;
1025     case 2:
1026     fRDCuts = dynamic_cast<AliRDHFCutsDStartoKpipi*> (GetOutputData(3));
1027     break;
1028     default:
1029     break;
1030     }
1031     if (!fRDCuts) {
1032     printf("ERROR: fRDCuts not available\n");
1033     return;
1034     }
1035
1036     // TCanvas* cvex=new TCanvas("example","Example Mass");
1037     // cvex->cd();
1038     // ((TH1F*)fOutput->FindObject("hMass_pt0phi0"))->Draw();
1039     // TCanvas* cv2d=new TCanvas("2d","mass-cos2phi");
1040     // cv2d->cd();
1041     // ((TH2F*)fOutput->FindObject("hMc2phi"))->Draw("colz");
1042     // TCanvas *cstat=new TCanvas("cstat","Stat");
1043     // cstat->SetGridy();
1044     // cstat->cd();
1045     // fhEventsInfo->Draw("htext0");
1046
1047     //  TMultiGraph *multig = new TMultiGraph();
1048     if(fDecChannel==2)return;
1049     TGraphErrors *g[fNPtBins];
1050     TH1F *h = new TH1F("h","h",100,0.,1.);
1051     TString hname;
1052     TString gname;
1053     for(Int_t ipt = 0;ipt<fNPtBins;ipt++){
1054     g[ipt] = new TGraphErrors(fNPhiBinLims);
1055     gname.Form("hMass_pt%d",ipt);
1056     g[ipt]->SetTitle(gname.Data());
1057     for(Int_t iphi = 0;iphi<fNPhiBinLims;iphi++){
1058     hname.Form("hMass_pt%dphi%d",ipt,iphi);
1059     h = (TH1F*)fOutput->FindObject("hMass_pt0phi0");
1060     AliHFMassFitter fitter(h,fLowmasslimit,fUpmasslimit,2,0);
1061     Int_t pdg=0;
1062     switch(fDecChannel){
1063     case 0:
1064     pdg=411;
1065     break;
1066     case 1:
1067     pdg=421;
1068     break;
1069     case 2:
1070     pdg=413;
1071     break;
1072     default:
1073     break;
1074     }
1075     fitter.SetInitialGaussianMean(TDatabasePDG::Instance()->GetParticle(pdg)->Mass());
1076     fitter.SetInitialGaussianSigma(0.012);
1077     fitter.InitNtuParam(Form("ntuPtbin%d",ipt));
1078     Double_t signal=0, errSignal=0;
1079     if(fitter.MassFitter(kFALSE)){
1080     fitter.Signal(3,signal,errSignal);
1081     }
1082     g[ipt]->SetPoint(iphi,fPhiBins[iphi],signal);
1083     g[ipt]->SetPointError(iphi,fPhiBins[iphi],errSignal);
1084     }//end loop over phi
1085      //     multig->Add(g[ipt],"ap");
1086      }//end loop on pt
1087      TCanvas *cdndphi = new TCanvas("dN/d#phi","dN/d#phi");
1088      cdndphi->Divide(1,fNPtBins);
1089      for(Int_t ipt = 0;ipt<fNPtBins;ipt++){
1090      cdndphi->cd(ipt+1);
1091      g[ipt]->Draw("AP");
1092      }
1093   */
1094   return;
1095 }
1096
1097 //-------------------------------------------
1098 /*
1099   Float_t GetEventPlaneFromVZERO(){
1100   AliAODHFUtil *info = (AliAODHFUtil*) aodEvent->GetList()->FindObject("fHFUtilInfo");
1101   if (!info) return -999.;
1102   //============= FIX ONLY FOR AOD033
1103   Double_t *par0;
1104   Double_t par0_137161[64] = { 6.71e-02 , 6.86e-02 , 7.06e-02 , 6.32e-02 ,
1105   5.91e-02 , 6.07e-02 , 5.78e-02 , 5.73e-02 , 5.91e-02 , 6.22e-02 ,
1106   5.90e-02 , 6.11e-02 , 5.55e-02 , 5.29e-02 , 5.19e-02 , 5.56e-02 ,
1107   6.25e-02 , 7.03e-02 , 5.64e-02 , 5.81e-02 , 4.57e-02 , 5.30e-02 ,
1108   5.13e-02 , 6.43e-02 , 6.27e-02 , 6.48e-02 , 6.07e-02 , 1.01e-01 ,
1109   6.68e-02 , 7.16e-02 , 6.36e-02 , 5.95e-02 , 2.52e-02 , 2.82e-02 ,
1110   2.56e-02 , 2.86e-02 , 2.82e-02 , 2.10e-02 , 2.13e-02 , 2.32e-02 ,
1111   2.75e-02 , 4.34e-02 , 3.78e-02 , 4.52e-02 , 4.11e-02 , 3.89e-02 ,
1112   4.10e-02 , 3.73e-02 , 4.51e-02 , 5.07e-02 , 5.42e-02 , 4.74e-02 ,
1113   4.33e-02 , 4.44e-02 , 4.64e-02 , 3.01e-02 , 6.38e-02 , 5.26e-02 ,
1114   4.99e-02 , 5.26e-02 , 5.47e-02 , 3.84e-02 , 5.00e-02 , 5.20e-02 };
1115   Double_t par0_137366[64] = { 7.12e-02 , 7.34e-02 , 7.39e-02 , 6.54e-02 , 6.11e-02 , 6.31e-02 , 6.15e-02 ,
1116   6.00e-02 , 6.10e-02 , 6.49e-02 , 6.17e-02 , 6.33e-02 , 6.00e-02 , 5.48e-02 ,
1117   5.44e-02 , 5.81e-02 , 6.49e-02 , 7.07e-02 , 5.91e-02 , 6.18e-02 , 4.82e-02 ,
1118   5.67e-02 , 5.36e-02 , 6.60e-02 , 6.37e-02 , 6.78e-02 , 6.31e-02 , 1.04e-01 ,
1119   6.91e-02 , 7.32e-02 , 6.61e-02 , 6.16e-02 , 2.64e-02 , 2.81e-02 , 2.64e-02 ,
1120   2.85e-02 , 2.87e-02 , 2.18e-02 , 2.19e-02 , 2.43e-02 , 2.81e-02 , 4.37e-02 ,
1121   3.90e-02 , 4.66e-02 , 4.24e-02 , 4.09e-02 , 4.21e-02 , 3.88e-02 , 4.83e-02 ,
1122   5.23e-02 , 5.44e-02 , 4.85e-02 , 4.42e-02 , 4.58e-02 , 4.74e-02 , 3.14e-02 ,
1123   6.31e-02 , 5.30e-02 , 5.01e-02 , 5.33e-02 , 5.70e-02 , 3.95e-02 , 4.98e-02 , 5.31e-02 };
1124   if (aodEvent->GetRunNumber() <= 137165)
1125   par0=par0_137161;
1126   else
1127   par0=par0_137366;
1128   Float_t vChCorr[64];
1129   for(int i=0; i!=64; ++i)
1130   vChCorr[i] = (info->GetVZEROChannel(i))/par0[i]/64.;
1131   //============= END OF FIX AOD033
1132   Float_t multR[8];
1133   for(int i=0; i!=8; ++i) multR[i] = 0;
1134   for(int i=0; i!=64; ++i)
1135   multR[i/8] += vChCorr[i];
1136   for(int i=0; i!=8; ++i)
1137   if(multR[i]) {
1138   double Qx=0, Qy=0;
1139   for(int j=0; j!=8; ++j) {
1140   double phi = TMath::PiOver4()*(0.5+j);
1141   Qx += TMath::Cos(2*phi)*vChCorr[i*8+j]/multR[i];
1142   Qy += TMath::Sin(2*phi)*vChCorr[i*8+j]/multR[i];
1143   }
1144   }
1145   return 0.5*TMath::ATan2(Qy,Qx)+TMath::PiOver2();
1146   }
1147
1148 */