[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 "AliAODEvent.h"
#include "AliAODVertex.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);

  fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec");

  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());
  
  fHistNEvents->Fill(0); // count event
  // Post the data already here
  PostData(1,fOutput);

  // load HF vertices                     
  TClonesArray *inputArrayVertices =
    (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
  if(!inputArrayVertices) {
    printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
    return;
  }

  // load D0->Kpi candidates                                                   
  TClonesArray *inputArrayD0toKpi =
    (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
  if(!inputArrayD0toKpi) {
    printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
    return;
  }

 
  // load D*+ candidates                                                   
  TClonesArray *inputArrayDstar =
    (TClonesArray*)aod->GetList()->FindObject("Dstar");
  if(!inputArrayDstar) {
    printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
    return;
  }
  

  // 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,cov[6],errx,erry,errz;
  AliAODRecoDecayHF2Prong *d2=0;
  AliAODRecoDecayHF3Prong *d3=0;
  AliAODRecoDecayHF4Prong *d4=0;

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

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

    // get parent AliAODRecoDecayHF
    nprongs= vtx->GetNDaughters();

    switch(nprongs) {
    case 2: // look for D0->Kpi
      d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
      if(d2->Charge() != 0) continue; // these are D* 
      lab = d2->MatchToMC(421,mcArray);
      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)) {
	  AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(lab);
	  pdg = part->GetPdgCode();
	  invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
	  fHistMass->Fill(invmass);
	  // Post the data already here
	  PostData(1,fOutput);
	  
	  AliAODVertex *secVtx = d2->GetSecondaryVtx();
	  secVtx->GetCovarianceMatrix(cov);
	  errx=1.; erry=1.; errz=1.;
	  if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
	  if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
	  if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
	  
	  fNtupleCmp->Fill(pdg,nprongs,d2->Xv(),part->Xv(),errx,d2->Yv(),part->Yv(),erry,d2->Zv(),part->Zv(),errz,d2->GetReducedChi2(),d2->Pt(),invmass);
	  PostData(2,fNtupleCmp);
	}
	if(unsetvtx) d2->UnsetOwnPrimaryVtx();
      }
      break;
    case 3: // look for D+
      d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
      lab = d3->MatchToMC(411,mcArray);
      if(lab>=0) {
	unsetvtx=kFALSE;
	if(!d3->GetOwnPrimaryVtx()) {
	  d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	  unsetvtx=kTRUE;
	}
	if(d3->SelectDplus(fVHF->GetDplusCuts())) {
	  AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(lab);
	  pdg = part->GetPdgCode();
	  invmass = d3->InvMassDplus();
	  AliAODVertex *secVtx = d3->GetSecondaryVtx();
	  secVtx->GetCovarianceMatrix(cov);
	  errx=1.; erry=1.; errz=1.;
	  if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
	  if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
	  if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
	  
	  fNtupleCmp->Fill(pdg,nprongs,d3->Xv(),part->Xv(),errx,d3->Yv(),part->Yv(),erry,d3->Zv(),part->Zv(),errz,d3->GetReducedChi2(),d3->Pt(),invmass);
	  PostData(2,fNtupleCmp);
	}
	if(unsetvtx) d3->UnsetOwnPrimaryVtx();
      }
      break;
    case 4: // look for D0->Kpipipi
      d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
      lab = d4->MatchToMC(421,mcArray);
      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 *part = (AliAODMCParticle*)mcArray->At(lab);
	  pdg = part->GetPdgCode();
	  //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
	  invmass = 10.; 	  

	  AliAODVertex *secVtx = d4->GetSecondaryVtx();
	  secVtx->GetCovarianceMatrix(cov);
	  errx=1.; erry=1.; errz=1.;
	  if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
	  if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
	  if(cov[5]>0) errz = TMath::Sqrt(cov[5]);

	  fNtupleCmp->Fill(pdg,nprongs,d4->Xv(),part->Xv(),errx,d4->Yv(),part->Yv(),erry,d4->Zv(),part->Zv(),errz,d4->GetReducedChi2(),d4->Pt(),invmass);
	  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
	29	#include "AliAODEvent.h"
	30	#include "AliAODVertex.h"
	31	#include "AliAODTrack.h"
	32	#include "AliAODMCHeader.h"
	33	#include "AliAODMCParticle.h"
	34	#include "AliAODRecoDecayHF2Prong.h"
b3999999	35	#include "AliAODRecoCascadeHF.h"
2ae9c7cb	36	#include "AliAnalysisVertexingHF.h"
3a219f60	37	#include "AliAnalysisTaskSE.h"
	38	#include "AliAnalysisTaskSECompareHF.h"
	39
	40	ClassImp(AliAnalysisTaskSECompareHF)
	41
	42
	43	//________________________________________________________________________
	44	AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
	45	AliAnalysisTaskSE(),
	46	fOutput(0),
2a041947	47	fNtupleCmp(0),
12bfc069	48	fHistMass(0),
a273622f	49	fHistNEvents(0),
12bfc069	50	fVHF(0)
3a219f60	51	{
3a219f60	52	// Default constructor
3a219f60	53
9d6d35b0	54	// NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
3a219f60	55	}
	56
	57	//________________________________________________________________________
	58	AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
	59	AliAnalysisTaskSE(name),
	60	fOutput(0),
2a041947	61	fNtupleCmp(0),
12bfc069	62	fHistMass(0),
a273622f	63	fHistNEvents(0),
12bfc069	64	fVHF(0)
3a219f60	65	{
9d6d35b0	66	// Standard constructor
3a219f60	67
	68	// Output slot #1 writes into a TList container
	69	DefineOutput(1,TList::Class()); //My private output
10bdd1ae	70	// Output slot #2 writes into a TNtuple container
10bdd1ae	71	DefineOutput(2,TNtuple::Class()); //My private output
3a219f60	72	}
	73
	74	//________________________________________________________________________
	75	AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
	76	{
	77	// Destructor
	78	if (fOutput) {
	79	delete fOutput;
	80	fOutput = 0;
	81	}
12bfc069	82	if (fVHF) {
	83	delete fVHF;
	84	fVHF = 0;
	85	}
3a219f60	86	}
	87
	88	//________________________________________________________________________
	89	void AliAnalysisTaskSECompareHF::Init()
	90	{
	91	// Initialization
	92
	93	if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
9d6d35b0	94
12bfc069	95	gROOT->LoadMacro("ConfigVertexingHF.C");
	96
	97	fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
	98	fVHF->PrintStatus();
9d6d35b0	99
3a219f60	100	return;
	101	}
	102
	103	//________________________________________________________________________
	104	void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
	105	{
	106	// Create the output container
	107	//
	108	if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
	109
	110	// Several histograms are more conveniently managed in a TList
	111	fOutput = new TList();
	112	fOutput->SetOwner();
	113
	114	fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
	115	fHistMass->Sumw2();
	116	fHistMass->SetMinimum(0);
	117	fOutput->Add(fHistMass);
	118
a273622f	119	fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
	120	fHistNEvents->Sumw2();
	121	fHistNEvents->SetMinimum(0);
	122	fOutput->Add(fHistNEvents);
	123
2a041947	124	fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec");
3a219f60	125
	126	return;
	127	}
	128
	129	//________________________________________________________________________
	130	void AliAnalysisTaskSECompareHF::UserExec(Option_t /option*/)
	131	{
	132	// Execute analysis for current event:
	133	// heavy flavor candidates association to MC truth
9d6d35b0	134
3a219f60	135
3a219f60	136	AliAODEvent aod = dynamic_cast<AliAODEvent> (InputEvent());
9d6d35b0	137
a273622f	138	fHistNEvents->Fill(0); // count event
	139	// Post the data already here
	140	PostData(1,fOutput);
2a041947	141
	142	// load HF vertices
	143	TClonesArray *inputArrayVertices =
	144	(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
	145	if(!inputArrayVertices) {
	146	printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
	147	return;
	148	}
	149
3a219f60	150	// load D0->Kpi candidates
	151	TClonesArray *inputArrayD0toKpi =
	152	(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
	153	if(!inputArrayD0toKpi) {
	154	printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
	155	return;
	156	}
	157
9d6d35b0	158
b3999999	159	// load D*+ candidates
	160	TClonesArray *inputArrayDstar =
	161	(TClonesArray*)aod->GetList()->FindObject("Dstar");
	162	if(!inputArrayDstar) {
	163	printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
	164	return;
	165	}
9d6d35b0	166
b3999999	167
3a219f60	168	// AOD primary vertex
	169	AliAODVertex vtx1 = (AliAODVertex)aod->GetPrimaryVertex();
	170	//vtx1->Print();
	171
	172	// load MC particles
	173	TClonesArray *mcArray =
	174	(TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
	175	if(!mcArray) {
	176	printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
	177	return;
	178	}
	179
	180	// load MC header
	181	AliAODMCHeader *mcHeader =
	182	(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
	183	if(!mcHeader) {
	184	printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
	185	return;
	186	}
	187
2a041947	188	Int_t nprongs,lab,okD0,okD0bar,pdg;
	189	Bool_t unsetvtx;
	190	Double_t invmass,cov[6],errx,erry,errz;
952f1e33	191	AliAODRecoDecayHF2Prong *d2=0;
	192	AliAODRecoDecayHF3Prong *d3=0;
	193	AliAODRecoDecayHF4Prong *d4=0;
2a041947	194
	195	// loop over vertices
	196	Int_t nVertices = inputArrayVertices->GetEntriesFast();
	197	printf("Number of vertices: %d\n",nVertices);
	198
	199	for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
	200	AliAODVertex vtx = (AliAODVertex)inputArrayVertices->UncheckedAt(iVtx);
	201
	202	// get parent AliAODRecoDecayHF
	203	nprongs= vtx->GetNDaughters();
	204
	205	switch(nprongs) {
	206	case 2: // look for D0->Kpi
952f1e33	207	d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
a273622f	208	if(d2->Charge() != 0) continue; // these are D*
2a041947	209	lab = d2->MatchToMC(421,mcArray);
	210	if(lab>=0) {
	211	unsetvtx=kFALSE;
	212	if(!d2->GetOwnPrimaryVtx()) {
	213	d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	214	unsetvtx=kTRUE;
	215	}
	216	okD0=0; okD0bar=0;
	217	if(d2->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {
	218	AliAODMCParticle part = (AliAODMCParticle)mcArray->At(lab);
	219	pdg = part->GetPdgCode();
	220	invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
	221	fHistMass->Fill(invmass);
	222	// Post the data already here
	223	PostData(1,fOutput);
	224
	225	AliAODVertex *secVtx = d2->GetSecondaryVtx();
	226	secVtx->GetCovarianceMatrix(cov);
	227	errx=1.; erry=1.; errz=1.;
	228	if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
	229	if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
	230	if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
	231
	232	fNtupleCmp->Fill(pdg,nprongs,d2->Xv(),part->Xv(),errx,d2->Yv(),part->Yv(),erry,d2->Zv(),part->Zv(),errz,d2->GetReducedChi2(),d2->Pt(),invmass);
	233	PostData(2,fNtupleCmp);
	234	}
	235	if(unsetvtx) d2->UnsetOwnPrimaryVtx();
3a219f60	236	}
2a041947	237	break;
2a041947	238	case 3: // look for D+
952f1e33	239	d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
2a041947	240	lab = d3->MatchToMC(411,mcArray);
	241	if(lab>=0) {
	242	unsetvtx=kFALSE;
	243	if(!d3->GetOwnPrimaryVtx()) {
	244	d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	245	unsetvtx=kTRUE;
	246	}
	247	if(d3->SelectDplus(fVHF->GetDplusCuts())) {
	248	AliAODMCParticle part = (AliAODMCParticle)mcArray->At(lab);
	249	pdg = part->GetPdgCode();
	250	invmass = d3->InvMassDplus();
	251	AliAODVertex *secVtx = d3->GetSecondaryVtx();
	252	secVtx->GetCovarianceMatrix(cov);
	253	errx=1.; erry=1.; errz=1.;
	254	if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
	255	if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
	256	if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
	257
	258	fNtupleCmp->Fill(pdg,nprongs,d3->Xv(),part->Xv(),errx,d3->Yv(),part->Yv(),erry,d3->Zv(),part->Zv(),errz,d3->GetReducedChi2(),d3->Pt(),invmass);
	259	PostData(2,fNtupleCmp);
	260	}
	261	if(unsetvtx) d3->UnsetOwnPrimaryVtx();
	262	}
	263	break;
	264	case 4: // look for D0->Kpipipi
952f1e33	265	d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
2a041947	266	lab = d4->MatchToMC(421,mcArray);
	267	if(lab>=0) {
	268	unsetvtx=kFALSE;
	269	if(!d4->GetOwnPrimaryVtx()) {
	270	d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	271	unsetvtx=kTRUE;
	272	}
	273	okD0=0; okD0bar=0;
	274	if(d4->SelectD0(fVHF->GetD0to4ProngsCuts(),okD0,okD0bar)) {
	275	AliAODMCParticle part = (AliAODMCParticle)mcArray->At(lab);
	276	pdg = part->GetPdgCode();
	277	//invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
	278	invmass = 10.;
	279
	280	AliAODVertex *secVtx = d4->GetSecondaryVtx();
	281	secVtx->GetCovarianceMatrix(cov);
	282	errx=1.; erry=1.; errz=1.;
	283	if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
	284	if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
	285	if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
	286
	287	fNtupleCmp->Fill(pdg,nprongs,d4->Xv(),part->Xv(),errx,d4->Yv(),part->Yv(),erry,d4->Zv(),part->Zv(),errz,d4->GetReducedChi2(),d4->Pt(),invmass);
	288	PostData(2,fNtupleCmp);
	289
	290	}
	291	if(unsetvtx) d4->UnsetOwnPrimaryVtx();
	292	}
	293	break;
3a219f60	294	}
3a219f60	295
2a041947	296	} // end loop on vertices
3a219f60	297
9d6d35b0	298
b3999999	299	// loop over D*+ candidates
4dd17b11	300	/*
b3999999	301	for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
	302	AliAODRecoCascadeHF c = (AliAODRecoCascadeHF)inputArrayDstar->UncheckedAt(iDstar);
	303	Int_t labDstar = c->MatchToMC(413,421,mcArray);
	304	if(labDstar>=0) printf("GOOD MATCH FOR D*+\n");
	305	}
4dd17b11	306	*/
2a041947	307
4dd17b11	308	// Post the data already here
4dd17b11	309	PostData(1,fOutput);
2a041947	310	PostData(2,fNtupleCmp);
b3999999	311
3a219f60	312	return;
3a219f60	313	}
3a219f60	314	//________________________________________________________________________
	315	void AliAnalysisTaskSECompareHF::Terminate(Option_t /option*/)
	316	{
	317	// Terminate analysis
	318	//
	319	if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
	320
	321	fOutput = dynamic_cast<TList*> (GetOutputData(1));
	322	if (!fOutput) {
	323	printf("ERROR: fOutput not available\n");
	324	return;
	325	}
	326
3a219f60	327	fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
a273622f	328	fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
3a219f60	329
2a041947	330	fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));
10bdd1ae	331
3a219f60	332	return;
	333	}
	334