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