[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSECompareHF.cxx

/**************************************************************************
 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/////////////////////////////////////////////////////////////
//
// AliAnalysisTaskSE for the comparison of heavy flavor
// decay candidates with the MC truth.
//
// Author: A.Dainese, andrea.dainese@lnl.infn.it
/////////////////////////////////////////////////////////////

#include <TClonesArray.h>
#include <TNtuple.h>
#include <TList.h>
#include <TH1F.h>

#include "AliAnalysisManager.h"
#include "AliAODHandler.h"
#include "AliAODEvent.h"
#include "AliAODVertex.h"
#include "AliESDtrack.h"
#include "AliAODTrack.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliAODRecoDecayHF2Prong.h"
#include "AliAODRecoCascadeHF.h"
#include "AliAnalysisVertexingHF.h"
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisTaskSECompareHF.h"

ClassImp(AliAnalysisTaskSECompareHF)


//________________________________________________________________________
AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
AliAnalysisTaskSE(),
fOutput(0), 
fNtupleCmp(0),
fHistMass(0),
fHistNEvents(0),
fVHF(0)
{
  // Default constructor

  // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
}

//________________________________________________________________________
AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
AliAnalysisTaskSE(name),
fOutput(0), 
fNtupleCmp(0),
fHistMass(0),
fHistNEvents(0),
fVHF(0)
{
  // Standard constructor

  // Output slot #1 writes into a TList container
  DefineOutput(1,TList::Class());  //My private output
  // Output slot #2 writes into a TNtuple container
  DefineOutput(2,TNtuple::Class());  //My private output
}

//________________________________________________________________________
AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
{
  // Destructor
  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
  if (fVHF) {
    delete fVHF;
    fVHF = 0;
  }
}  

//________________________________________________________________________
void AliAnalysisTaskSECompareHF::Init()
{
  // Initialization

  if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
  
  gROOT->LoadMacro("ConfigVertexingHF.C");

  fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  
  fVHF->PrintStatus();
  
  return;
}

//________________________________________________________________________
void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
{
  // Create the output container
  //
  if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");

  // Several histograms are more conveniently managed in a TList
  fOutput = new TList();
  fOutput->SetOwner();

  fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
  fHistMass->Sumw2();
  fHistMass->SetMinimum(0);
  fOutput->Add(fHistMass);

  fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
  fHistNEvents->Sumw2();
  fHistNEvents->SetMinimum(0);
  fOutput->Add(fHistNEvents);

  OpenFile(2); // 2 is the slot number of the ntuple
  fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec:CPta:Prodd0");

  return;
}

//________________________________________________________________________
void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
{
  // Execute analysis for current event:
  // heavy flavor candidates association to MC truth

  
  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());

  TClonesArray *inputArrayVertices = 0;
  TClonesArray *inputArrayD0toKpi = 0;
  TClonesArray *inputArrayDstar = 0;

  if(!aod && AODEvent() && IsStandardAOD()) {
    // In case there is an AOD handler writing a standard AOD, use the AOD 
    // event in memory rather than the input (ESD) event.    
    aod = dynamic_cast<AliAODEvent*> (AODEvent());
    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
    // have to taken from the AOD event hold by the AliAODExtension
    AliAODHandler* aodHandler = (AliAODHandler*) 
      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
    if(aodHandler->GetExtensions()) {
      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
      AliAODEvent *aodFromExt = ext->GetAOD();
      // load HF vertices                
      inputArrayVertices = (TClonesArray*)aodFromExt->GetList()->FindObject("VerticesHF");
      // load D0->Kpi candidates
      inputArrayD0toKpi = (TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
      // load D*+ candidates                                                   
      inputArrayDstar = (TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
    }
  } else {
    // load HF vertices                
    inputArrayVertices = (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
    // load D0->Kpi candidates                                                 
    inputArrayD0toKpi = (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
    // load D*+ candidates                                                   
    inputArrayDstar = (TClonesArray*)aod->GetList()->FindObject("Dstar");
  }


  if(!inputArrayVertices) {
    printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
    return;
  }
  if(!inputArrayD0toKpi) {
    printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
    return;
  }
  if(!inputArrayDstar) {
    printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
    return;
  }
  

  fHistNEvents->Fill(0); // count event
  // Post the data already here
  PostData(1,fOutput);

  // AOD primary vertex
  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
  //vtx1->Print();

  // load MC particles
  TClonesArray *mcArray = 
    (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
  if(!mcArray) {
    printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
    return;
  }

  // load MC header
  AliAODMCHeader *mcHeader = 
    (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
  if(!mcHeader) {
    printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
    return;
  }
  
  Int_t nprongs,lab,okD0,okD0bar,pdg;
  Bool_t unsetvtx;    
  Double_t invmass,posRec[3],posTrue[3],covRec[6],errx,erry,errz;
  AliAODRecoDecayHF2Prong *d2=0;
  AliAODRecoDecayHF3Prong *d3=0;
  AliAODRecoDecayHF4Prong *d4=0;

  Int_t pdgDgD0toKpi[2]={321,211};
  Int_t pdgDgDplustoKpipi[3]={321,211,211};
  Int_t pdgDgD0toKpipipi[4]={321,211,211,211};

  // loop over vertices
  Int_t nVertices = inputArrayVertices->GetEntriesFast();
  if(fDebug>1) printf("Number of vertices: %d\n",nVertices);

  for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
    AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx);

    vtx->GetXYZ(posRec);  
    vtx->GetCovarianceMatrix(covRec);
    errx=1.; erry=1.; errz=1.;
    if(covRec[0]>0) errx = TMath::Sqrt(covRec[0]);
    if(covRec[2]>0) erry = TMath::Sqrt(covRec[2]);
    if(covRec[5]>0) errz = TMath::Sqrt(covRec[5]);
	  

    // get parent AliAODRecoDecayHF
    nprongs= vtx->GetNDaughters();
    // check that the daughters have kITSrefit and kTPCrefit
    Bool_t allDgOK=kTRUE;
    for(Int_t idg=0; idg<nprongs; idg++) {
      AliAODTrack *track = (AliAODTrack*)vtx->GetDaughter(idg);
      if(!(track->GetStatus()&AliESDtrack::kITSrefit)) allDgOK=kFALSE;
      if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) allDgOK=kFALSE;
    }
    if(!allDgOK) continue;


    switch(nprongs) {
    case 2: // look for D0->Kpi
      d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
      if(d2->IsLikeSign()) continue;
      if(d2->Charge() != 0) continue; // these are D* 
      lab = d2->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
      if(lab>=0) {
	unsetvtx=kFALSE;
	if(!d2->GetOwnPrimaryVtx()) {
	  d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	  unsetvtx=kTRUE;
	}
	okD0=0; okD0bar=0; 
	if(d2->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) { // beware: cuts may bias the resolution!
	  AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
	  pdg = dMC->GetPdgCode();
	  invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
	  // get a daughter for true pos of decay vertex
	  AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
	  dg0MC->XvYvZv(posTrue);
	  fHistMass->Fill(invmass);
	  // Post the data already here
	  PostData(1,fOutput);
	  Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
			   (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
			   (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
			   (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
			   (Float_t)d2->GetReducedChi2(),(Float_t)d2->Pt(),(Float_t)invmass,
			   (Float_t)d2->CosPointingAngle(),(Float_t)d2->Prodd0d0()};
	  fNtupleCmp->Fill(tmp);
	  PostData(2,fNtupleCmp);
	}
	if(unsetvtx) d2->UnsetOwnPrimaryVtx();
      }
      break;
    case 3: // look for D+
      d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
      if(d3->IsLikeSign()) continue;
      lab = d3->MatchToMC(411,mcArray,3,pdgDgDplustoKpipi);
      if(lab>=0) {
	unsetvtx=kFALSE;
	if(!d3->GetOwnPrimaryVtx()) {
	  d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	  unsetvtx=kTRUE;
	}
	//if(d3->SelectDplus(fVHF->GetDplusCuts())) {
	AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
	pdg = dMC->GetPdgCode();
	invmass = d3->InvMassDplus();
	// get a daughter for true pos of decay vertex
	AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
	dg0MC->XvYvZv(posTrue);
	Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
			 (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
			 (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
			 (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
			 (Float_t)d3->GetReducedChi2(),(Float_t)d3->Pt(),(Float_t)invmass,
			 (Float_t)d3->CosPointingAngle(),(Float_t)(d3->Getd0Prong(0)*d3->Getd0Prong(1)*d3->Getd0Prong(2))};
	fNtupleCmp->Fill(tmp);
	PostData(2,fNtupleCmp);
	//}
	if(unsetvtx) d3->UnsetOwnPrimaryVtx();
      }
      break;
    case 4: // look for D0->Kpipipi
      d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
      if(d4->IsLikeSign()) continue;
      lab = d4->MatchToMC(421,mcArray,4,pdgDgD0toKpipipi);
      if(lab>=0) {
	unsetvtx=kFALSE;
	if(!d4->GetOwnPrimaryVtx()) {
	  d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	  unsetvtx=kTRUE;
	}
	okD0=0; okD0bar=0; 
	//if(d4->SelectD0(fVHF->GetD0to4ProngsCuts(),okD0,okD0bar)) {
	AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
	pdg = dMC->GetPdgCode();
	//invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
	invmass = 10.; 	  // dummy
	// get a daughter for true pos of decay vertex
	AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
	dg0MC->XvYvZv(posTrue);
	Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
			 (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
			 (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
			 (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
			 (Float_t)d4->GetReducedChi2(),(Float_t)d4->Pt(),(Float_t)invmass,
			 (Float_t)d4->CosPointingAngle(),(Float_t)(d4->Getd0Prong(0)*d4->Getd0Prong(1)*d4->Getd0Prong(2)*d4->Getd0Prong(3))};
	fNtupleCmp->Fill(tmp);
	PostData(2,fNtupleCmp);
	//}
	if(unsetvtx) d4->UnsetOwnPrimaryVtx();
      }
      break;
    }

  } // end loop on vertices


  // loop over D*+ candidates
  /*
  for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
    AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
    Int_t labDstar = c->MatchToMC(413,421,mcArray);
    if(labDstar>=0) printf("GOOD MATCH FOR D*+\n"); 
  }
  */

  // Post the data already here
  PostData(1,fOutput);
  PostData(2,fNtupleCmp);

  return;
}
//________________________________________________________________________
void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  //
  if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");

  fOutput = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutput) {     
    printf("ERROR: fOutput not available\n");
    return;
  }

  fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));

  //fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));

  return;
}
Commit	Line	Data
3a219f60	1	/**************************************************************************
	2	* Copyright(c) 1998-2008, 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	// AliAnalysisTaskSE for the comparison of heavy flavor
	19	// decay candidates with the MC truth.
	20	//
	21	// Author: A.Dainese, andrea.dainese@lnl.infn.it
	22	/////////////////////////////////////////////////////////////
	23
	24	#include <TClonesArray.h>
	25	#include <TNtuple.h>
	26	#include <TList.h>
	27	#include <TH1F.h>
	28
b557eb43	29	#include "AliAnalysisManager.h"
b557eb43	30	#include "AliAODHandler.h"
3a219f60	31	#include "AliAODEvent.h"
3a219f60	32	#include "AliAODVertex.h"
607719af	33	#include "AliESDtrack.h"
3a219f60	34	#include "AliAODTrack.h"
	35	#include "AliAODMCHeader.h"
	36	#include "AliAODMCParticle.h"
	37	#include "AliAODRecoDecayHF2Prong.h"
b3999999	38	#include "AliAODRecoCascadeHF.h"
2ae9c7cb	39	#include "AliAnalysisVertexingHF.h"
3a219f60	40	#include "AliAnalysisTaskSE.h"
	41	#include "AliAnalysisTaskSECompareHF.h"
	42
	43	ClassImp(AliAnalysisTaskSECompareHF)
	44
	45
	46	//________________________________________________________________________
	47	AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
	48	AliAnalysisTaskSE(),
	49	fOutput(0),
2a041947	50	fNtupleCmp(0),
12bfc069	51	fHistMass(0),
a273622f	52	fHistNEvents(0),
12bfc069	53	fVHF(0)
3a219f60	54	{
3a219f60	55	// Default constructor
3a219f60	56
9d6d35b0	57	// NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
3a219f60	58	}
	59
	60	//________________________________________________________________________
	61	AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
	62	AliAnalysisTaskSE(name),
	63	fOutput(0),
2a041947	64	fNtupleCmp(0),
12bfc069	65	fHistMass(0),
a273622f	66	fHistNEvents(0),
12bfc069	67	fVHF(0)
3a219f60	68	{
9d6d35b0	69	// Standard constructor
3a219f60	70
	71	// Output slot #1 writes into a TList container
	72	DefineOutput(1,TList::Class()); //My private output
10bdd1ae	73	// Output slot #2 writes into a TNtuple container
10bdd1ae	74	DefineOutput(2,TNtuple::Class()); //My private output
3a219f60	75	}
	76
	77	//________________________________________________________________________
	78	AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
	79	{
	80	// Destructor
	81	if (fOutput) {
	82	delete fOutput;
	83	fOutput = 0;
	84	}
12bfc069	85	if (fVHF) {
	86	delete fVHF;
	87	fVHF = 0;
	88	}
3a219f60	89	}
	90
	91	//________________________________________________________________________
	92	void AliAnalysisTaskSECompareHF::Init()
	93	{
	94	// Initialization
	95
	96	if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
9d6d35b0	97
12bfc069	98	gROOT->LoadMacro("ConfigVertexingHF.C");
	99
	100	fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
	101	fVHF->PrintStatus();
9d6d35b0	102
3a219f60	103	return;
	104	}
	105
	106	//________________________________________________________________________
	107	void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
	108	{
	109	// Create the output container
	110	//
	111	if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
	112
	113	// Several histograms are more conveniently managed in a TList
	114	fOutput = new TList();
	115	fOutput->SetOwner();
	116
	117	fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
	118	fHistMass->Sumw2();
	119	fHistMass->SetMinimum(0);
	120	fOutput->Add(fHistMass);
	121
a273622f	122	fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
	123	fHistNEvents->Sumw2();
	124	fHistNEvents->SetMinimum(0);
	125	fOutput->Add(fHistNEvents);
	126
2a5c77c6	127	OpenFile(2); // 2 is the slot number of the ntuple
607719af	128	fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec:CPta:Prodd0");
3a219f60	129
	130	return;
	131	}
	132
	133	//________________________________________________________________________
	134	void AliAnalysisTaskSECompareHF::UserExec(Option_t /option*/)
	135	{
	136	// Execute analysis for current event:
	137	// heavy flavor candidates association to MC truth
9d6d35b0	138
3a219f60	139
3a219f60	140	AliAODEvent aod = dynamic_cast<AliAODEvent> (InputEvent());
8931c313	141
b557eb43	142	TClonesArray *inputArrayVertices = 0;
	143	TClonesArray *inputArrayD0toKpi = 0;
	144	TClonesArray *inputArrayDstar = 0;
	145
	146	if(!aod && AODEvent() && IsStandardAOD()) {
	147	// In case there is an AOD handler writing a standard AOD, use the AOD
	148	// event in memory rather than the input (ESD) event.
	149	aod = dynamic_cast<AliAODEvent*> (AODEvent());
	150	// in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
	151	// have to taken from the AOD event hold by the AliAODExtension
	152	AliAODHandler* aodHandler = (AliAODHandler*)
	153	((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
	154	if(aodHandler->GetExtensions()) {
	155	AliAODExtension ext = (AliAODExtension)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
	156	AliAODEvent *aodFromExt = ext->GetAOD();
	157	// load HF vertices
	158	inputArrayVertices = (TClonesArray*)aodFromExt->GetList()->FindObject("VerticesHF");
	159	// load D0->Kpi candidates
	160	inputArrayD0toKpi = (TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
	161	// load D*+ candidates
	162	inputArrayDstar = (TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
	163	}
	164	} else {
	165	// load HF vertices
	166	inputArrayVertices = (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
	167	// load D0->Kpi candidates
	168	inputArrayD0toKpi = (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
	169	// load D*+ candidates
	170	inputArrayDstar = (TClonesArray*)aod->GetList()->FindObject("Dstar");
	171	}
	172
2a041947	173
2a041947	174	if(!inputArrayVertices) {
	175	printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
	176	return;
	177	}
3a219f60	178	if(!inputArrayD0toKpi) {
	179	printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
	180	return;
	181	}
b3999999	182	if(!inputArrayDstar) {
	183	printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
	184	return;
	185	}
9d6d35b0	186
b3999999	187
b557eb43	188	fHistNEvents->Fill(0); // count event
	189	// Post the data already here
	190	PostData(1,fOutput);
	191
3a219f60	192	// AOD primary vertex
	193	AliAODVertex vtx1 = (AliAODVertex)aod->GetPrimaryVertex();
	194	//vtx1->Print();
	195
	196	// load MC particles
	197	TClonesArray *mcArray =
	198	(TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
	199	if(!mcArray) {
	200	printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
	201	return;
	202	}
	203
	204	// load MC header
	205	AliAODMCHeader *mcHeader =
	206	(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
	207	if(!mcHeader) {
	208	printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
	209	return;
	210	}
	211
2a041947	212	Int_t nprongs,lab,okD0,okD0bar,pdg;
2a041947	213	Bool_t unsetvtx;
90426506	214	Double_t invmass,posRec[3],posTrue[3],covRec[6],errx,erry,errz;
952f1e33	215	AliAODRecoDecayHF2Prong *d2=0;
	216	AliAODRecoDecayHF3Prong *d3=0;
	217	AliAODRecoDecayHF4Prong *d4=0;
2a041947	218
c8112d2f	219	Int_t pdgDgD0toKpi[2]={321,211};
	220	Int_t pdgDgDplustoKpipi[3]={321,211,211};
	221	Int_t pdgDgD0toKpipipi[4]={321,211,211,211};
	222
2a041947	223	// loop over vertices
2a041947	224	Int_t nVertices = inputArrayVertices->GetEntriesFast();
3b1598da	225	if(fDebug>1) printf("Number of vertices: %d\n",nVertices);
2a041947	226
	227	for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
	228	AliAODVertex vtx = (AliAODVertex)inputArrayVertices->UncheckedAt(iVtx);
	229
90426506	230	vtx->GetXYZ(posRec);
	231	vtx->GetCovarianceMatrix(covRec);
	232	errx=1.; erry=1.; errz=1.;
	233	if(covRec[0]>0) errx = TMath::Sqrt(covRec[0]);
	234	if(covRec[2]>0) erry = TMath::Sqrt(covRec[2]);
	235	if(covRec[5]>0) errz = TMath::Sqrt(covRec[5]);
	236
	237
2a041947	238	// get parent AliAODRecoDecayHF
2a041947	239	nprongs= vtx->GetNDaughters();
607719af	240	// check that the daughters have kITSrefit and kTPCrefit
	241	Bool_t allDgOK=kTRUE;
	242	for(Int_t idg=0; idg<nprongs; idg++) {
	243	AliAODTrack track = (AliAODTrack)vtx->GetDaughter(idg);
	244	if(!(track->GetStatus()&AliESDtrack::kITSrefit)) allDgOK=kFALSE;
	245	if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) allDgOK=kFALSE;
	246	}
	247	if(!allDgOK) continue;
	248
2a041947	249
	250	switch(nprongs) {
	251	case 2: // look for D0->Kpi
952f1e33	252	d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
27917274	253	if(d2->IsLikeSign()) continue;
a273622f	254	if(d2->Charge() != 0) continue; // these are D*
c8112d2f	255	lab = d2->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
2a041947	256	if(lab>=0) {
	257	unsetvtx=kFALSE;
	258	if(!d2->GetOwnPrimaryVtx()) {
	259	d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	260	unsetvtx=kTRUE;
	261	}
	262	okD0=0; okD0bar=0;
607719af	263	if(d2->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) { // beware: cuts may bias the resolution!
	264	AliAODMCParticle dMC = (AliAODMCParticle)mcArray->At(lab);
	265	pdg = dMC->GetPdgCode();
	266	invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
607719af	267	// get a daughter for true pos of decay vertex
	268	AliAODMCParticle dg0MC = (AliAODMCParticle)mcArray->At(dMC->GetDaughter(0));
	269	dg0MC->XvYvZv(posTrue);
caf4ca8b	270	fHistMass->Fill(invmass);
	271	// Post the data already here
	272	PostData(1,fOutput);
607719af	273	Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
	274	(Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
	275	(Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
	276	(Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
	277	(Float_t)d2->GetReducedChi2(),(Float_t)d2->Pt(),(Float_t)invmass,
	278	(Float_t)d2->CosPointingAngle(),(Float_t)d2->Prodd0d0()};
	279	fNtupleCmp->Fill(tmp);
	280	PostData(2,fNtupleCmp);
	281	}
2a041947	282	if(unsetvtx) d2->UnsetOwnPrimaryVtx();
3a219f60	283	}
2a041947	284	break;
2a041947	285	case 3: // look for D+
952f1e33	286	d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
27917274	287	if(d3->IsLikeSign()) continue;
c8112d2f	288	lab = d3->MatchToMC(411,mcArray,3,pdgDgDplustoKpipi);
2a041947	289	if(lab>=0) {
	290	unsetvtx=kFALSE;
	291	if(!d3->GetOwnPrimaryVtx()) {
	292	d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	293	unsetvtx=kTRUE;
	294	}
90426506	295	//if(d3->SelectDplus(fVHF->GetDplusCuts())) {
	296	AliAODMCParticle dMC = (AliAODMCParticle)mcArray->At(lab);
	297	pdg = dMC->GetPdgCode();
	298	invmass = d3->InvMassDplus();
	299	// get a daughter for true pos of decay vertex
	300	AliAODMCParticle dg0MC = (AliAODMCParticle)mcArray->At(dMC->GetDaughter(0));
	301	dg0MC->XvYvZv(posTrue);
607719af	302	Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
	303	(Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
	304	(Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
	305	(Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
	306	(Float_t)d3->GetReducedChi2(),(Float_t)d3->Pt(),(Float_t)invmass,
	307	(Float_t)d3->CosPointingAngle(),(Float_t)(d3->Getd0Prong(0)d3->Getd0Prong(1)d3->Getd0Prong(2))};
	308	fNtupleCmp->Fill(tmp);
a7c1934a	309	PostData(2,fNtupleCmp);
90426506	310	//}
2a041947	311	if(unsetvtx) d3->UnsetOwnPrimaryVtx();
	312	}
	313	break;
	314	case 4: // look for D0->Kpipipi
952f1e33	315	d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
27917274	316	if(d4->IsLikeSign()) continue;
c8112d2f	317	lab = d4->MatchToMC(421,mcArray,4,pdgDgD0toKpipipi);
2a041947	318	if(lab>=0) {
	319	unsetvtx=kFALSE;
	320	if(!d4->GetOwnPrimaryVtx()) {
	321	d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	322	unsetvtx=kTRUE;
	323	}
	324	okD0=0; okD0bar=0;
90426506	325	//if(d4->SelectD0(fVHF->GetD0to4ProngsCuts(),okD0,okD0bar)) {
	326	AliAODMCParticle dMC = (AliAODMCParticle)mcArray->At(lab);
	327	pdg = dMC->GetPdgCode();
	328	//invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
	329	invmass = 10.; // dummy
	330	// get a daughter for true pos of decay vertex
	331	AliAODMCParticle dg0MC = (AliAODMCParticle)mcArray->At(dMC->GetDaughter(0));
	332	dg0MC->XvYvZv(posTrue);
607719af	333	Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
	334	(Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
	335	(Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
	336	(Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
	337	(Float_t)d4->GetReducedChi2(),(Float_t)d4->Pt(),(Float_t)invmass,
	338	(Float_t)d4->CosPointingAngle(),(Float_t)(d4->Getd0Prong(0)d4->Getd0Prong(1)d4->Getd0Prong(2)*d4->Getd0Prong(3))};
	339	fNtupleCmp->Fill(tmp);
a7c1934a	340	PostData(2,fNtupleCmp);
90426506	341	//}
2a041947	342	if(unsetvtx) d4->UnsetOwnPrimaryVtx();
	343	}
	344	break;
3a219f60	345	}
3a219f60	346
2a041947	347	} // end loop on vertices
3a219f60	348
9d6d35b0	349
b3999999	350	// loop over D*+ candidates
4dd17b11	351	/*
b3999999	352	for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
	353	AliAODRecoCascadeHF c = (AliAODRecoCascadeHF)inputArrayDstar->UncheckedAt(iDstar);
	354	Int_t labDstar = c->MatchToMC(413,421,mcArray);
	355	if(labDstar>=0) printf("GOOD MATCH FOR D*+\n");
	356	}
4dd17b11	357	*/
2a041947	358
4dd17b11	359	// Post the data already here
4dd17b11	360	PostData(1,fOutput);
2a041947	361	PostData(2,fNtupleCmp);
b3999999	362
3a219f60	363	return;
3a219f60	364	}
3a219f60	365	//________________________________________________________________________
	366	void AliAnalysisTaskSECompareHF::Terminate(Option_t /option*/)
	367	{
	368	// Terminate analysis
	369	//
	370	if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
	371
	372	fOutput = dynamic_cast<TList*> (GetOutputData(1));
	373	if (!fOutput) {
	374	printf("ERROR: fOutput not available\n");
	375	return;
	376	}
	377
3a219f60	378	fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
a273622f	379	fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
3a219f60	380
27917274	381	//fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));
10bdd1ae	382
3a219f60	383	return;
	384	}
	385