[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDplus.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 extraction of signal(e.g D+) of heavy flavor
// decay candidates with the MC truth.
// Author: Renu Bala, bala@to.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 "AliAODRecoDecayHF3Prong.h"
#include "AliAnalysisVertexingHF.h"
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisTaskSEDplus.h"

ClassImp(AliAnalysisTaskSEDplus)


//________________________________________________________________________
AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
AliAnalysisTaskSE(),
fOutput(0), 
fNtupleDplus(0),
fNtupleDplusbackg(0),
fHistMass(0),
fHistSignal(0),
fHistBackground(0),
fVHF(0)
{
  // Default constructor
}

//________________________________________________________________________
AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name):
AliAnalysisTaskSE(name),
fOutput(0), 
fNtupleDplus(0),
fNtupleDplusbackg(0),
fHistMass(0),
fHistSignal(0),
fHistBackground(0),
fVHF(0)
{
  // Default 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
  // Output slot #3 writes into a TNtuple container
  DefineOutput(3,TNtuple::Class());  //My private output
}

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

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

  if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");

  gROOT->LoadMacro("ConfigVertexingHF.C");

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

  return;
}

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

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

  fHistMass = new TH1F("fHistMass", "D^{+} invariant mass; M [GeV]; Entries",100,1.765,1.965);
  fHistSignal = new TH1F("fHistSignal", "D^{+} invariant mass - MC; M [GeV]; Entries",100,1.765,1.965);
  fHistBackground =new TH1F("fHistBackground", "Background invariant mass - MC; M [GeV]; Entries",100,1.765,1.965);


  fHistMass->Sumw2();
  fHistSignal->Sumw2();
  fHistBackground->Sumw2();


  //fHistMass->SetMinimum(0);
  //fHistSignal->SetMinimum(0);
  //fHistBackground->SetMinimum(0);
  fOutput->Add(fHistMass);
  fOutput->Add(fHistSignal);
  fOutput->Add(fHistBackground);
  

  fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:Ptpi:PtK:Ptpi2:PtRec:PtTrue:PointingAngle:DecLeng:VxTrue:VxRec:VyRec,VzRec,InvMass:sigvert:d0Pi:d0K:d0Pi2");
  fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptkbkg:Ptpi2bkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:VyRecbkg:VzRecbkg:InvMassbkg:sigvertbkg:d0Pibkg:d0Kbkg:d0Pi2bkg");
  
  return;
}

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

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

  // load Dplus->Kpipi candidates                                                   
  TClonesArray *array3Prong =
    (TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
  if(!array3Prong) {
    printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
    return;
  }

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

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

  // load MC header
  AliAODMCHeader *mcHeader = 
    (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
  if(!mcHeader) {
    printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
    return;
  }
    Int_t n3Prong = array3Prong->GetEntriesFast();
    printf("Number of D+->Kpipi: %d\n",n3Prong);


    for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
      AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);


      Bool_t unsetvtx=kFALSE;
      if(!d->GetOwnPrimaryVtx()){
	d->SetOwnPrimaryVtx(vtx1);
	unsetvtx=kTRUE;
      }
      if(d->SelectDplus(fVHF->GetDplusCuts())) {
	  
	  
	Int_t labDp = d->MatchToMC(411,arrayMC);
    
	if(labDp>=0) {
	  AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
	  Int_t pdgDp = TMath::Abs(partDp->GetPdgCode());
	  if(pdgDp==411){
	    
	    fHistSignal->Fill(d->InvMassDplus());
	    fHistMass->Fill(d->InvMassDplus());
	    
	    // Post the data already here
	    PostData(1,fOutput);
	    
	    AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));

	    Float_t tmp[20];
	    tmp[0]=pdgDp;
	    tmp[1]=partDp->Px()-d->Px();
	    tmp[2]=partDp->Py()-d->Py();
	    tmp[3]=partDp->Pz()-d->Pz();
	    tmp[4]=d->PtProng(0);
	    tmp[5]=d->PtProng(1);
	    tmp[6]=d->PtProng(2);
	    tmp[7]=d->Pt();
	    tmp[8]=partDp->Pt();
	    tmp[9]=d->CosPointingAngle();
	    tmp[10]=d->DecayLength();
	    tmp[11]=dg0->Xv();
	    tmp[12]=d->Xv();
	    tmp[13]=d->Yv();
	    tmp[14]=d->Zv();
	    tmp[15]=d->InvMassDplus();
	    tmp[16]=d->GetSigmaVert();
	    tmp[17]=d->Getd0Prong(0);
	    tmp[18]=d->Getd0Prong(1);
	    tmp[19]=d->Getd0Prong(2);


	    fNtupleDplus->Fill(tmp);
	    PostData(2,fNtupleDplus);
	  }
	} else {     

	    Float_t tmpbkg[14];
	    tmpbkg[0]=d->PtProng(0);
	    tmpbkg[1]=d->PtProng(1);
	    tmpbkg[2]=d->PtProng(2);
	    tmpbkg[3]=d->Pt();
	    tmpbkg[4]=d->CosPointingAngle();
	    tmpbkg[5]=d->DecayLength();
	    tmpbkg[6]=d->Xv();
	    tmpbkg[7]=d->Yv();
	    tmpbkg[8]=d->Zv();
	    tmpbkg[9]=d->InvMassDplus();
	    tmpbkg[10]=d->GetSigmaVert();
	    tmpbkg[11]=d->Getd0Prong(0);
	    tmpbkg[12]=d->Getd0Prong(1);
	    tmpbkg[13]=d->Getd0Prong(2);

	  fHistBackground->Fill(d->InvMassDplus());
	  fNtupleDplusbackg->Fill(tmpbkg);	    
	  PostData(3,fNtupleDplusbackg);
	  fHistMass->Fill(d->InvMassDplus());
	}
	
	
      }

      if(unsetvtx) d->UnsetOwnPrimaryVtx();

    }
     // end loop on D+->Kpipi


  return;
}

//________________________________________________________________________
void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  //
  if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");

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

  fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
  fHistSignal = dynamic_cast<TH1F*>(fOutput->FindObject("fHistSignal"));
  fHistBackground = dynamic_cast<TH1F*>(fOutput->FindObject("fHistBackground"));
  fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
  fNtupleDplusbackg = dynamic_cast<TNtuple*>(GetOutputData(3));

  return;
}
Commit	Line	Data
d486095a	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 extraction of signal(e.g D+) of heavy flavor
	19	// decay candidates with the MC truth.
8c34bd86	20	// Author: Renu Bala, bala@to.infn.it
d486095a	21	/////////////////////////////////////////////////////////////
	22
	23	#include <TClonesArray.h>
	24	#include <TNtuple.h>
	25	#include <TList.h>
	26	#include <TH1F.h>
	27	#include "AliAODEvent.h"
	28	#include "AliAODVertex.h"
	29	#include "AliAODTrack.h"
	30	#include "AliAODMCHeader.h"
	31	#include "AliAODMCParticle.h"
	32	#include "AliAODRecoDecayHF3Prong.h"
	33	#include "AliAnalysisVertexingHF.h"
	34	#include "AliAnalysisTaskSE.h"
	35	#include "AliAnalysisTaskSEDplus.h"
	36
	37	ClassImp(AliAnalysisTaskSEDplus)
	38
	39
	40	//________________________________________________________________________
	41	AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
	42	AliAnalysisTaskSE(),
	43	fOutput(0),
	44	fNtupleDplus(0),
	45	fNtupleDplusbackg(0),
	46	fHistMass(0),
	47	fHistSignal(0),
	48	fHistBackground(0),
	49	fVHF(0)
	50	{
	51	// Default constructor
d486095a	52	}
	53
	54	//________________________________________________________________________
	55	AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name):
	56	AliAnalysisTaskSE(name),
	57	fOutput(0),
	58	fNtupleDplus(0),
	59	fNtupleDplusbackg(0),
	60	fHistMass(0),
	61	fHistSignal(0),
	62	fHistBackground(0),
	63	fVHF(0)
	64	{
	65	// Default constructor
	66
	67	// Output slot #1 writes into a TList container
	68	DefineOutput(1,TList::Class()); //My private output
10bdd1ae	69	// Output slot #2 writes into a TNtuple container
8c34bd86	70	DefineOutput(2,TNtuple::Class()); //My private output
10bdd1ae	71	// Output slot #3 writes into a TNtuple container
8c34bd86	72	DefineOutput(3,TNtuple::Class()); //My private output
d486095a	73	}
	74
	75	//________________________________________________________________________
	76	AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
	77	{
	78	// Destructor
	79	if (fOutput) {
	80	delete fOutput;
	81	fOutput = 0;
	82	}
	83	if (fVHF) {
	84	delete fVHF;
	85	fVHF = 0;
	86	}
	87	}
	88
	89	//________________________________________________________________________
	90	void AliAnalysisTaskSEDplus::Init()
	91	{
	92	// Initialization
	93
	94	if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
	95
	96	gROOT->LoadMacro("ConfigVertexingHF.C");
	97
	98	fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
	99	fVHF->PrintStatus();
	100
	101	return;
	102	}
	103
	104	//________________________________________________________________________
	105	void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
	106	{
	107	// Create the output container
	108	//
	109	if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
	110
	111	// Several histograms are more conveniently managed in a TList
	112	fOutput = new TList();
	113	fOutput->SetOwner();
	114
	115	fHistMass = new TH1F("fHistMass", "D^{+} invariant mass; M [GeV]; Entries",100,1.765,1.965);
	116	fHistSignal = new TH1F("fHistSignal", "D^{+} invariant mass - MC; M [GeV]; Entries",100,1.765,1.965);
	117	fHistBackground =new TH1F("fHistBackground", "Background invariant mass - MC; M [GeV]; Entries",100,1.765,1.965);
	118
	119
	120
	121
	122	fHistMass->Sumw2();
	123	fHistSignal->Sumw2();
	124	fHistBackground->Sumw2();
	125
	126
	127	//fHistMass->SetMinimum(0);
	128	//fHistSignal->SetMinimum(0);
	129	//fHistBackground->SetMinimum(0);
	130	fOutput->Add(fHistMass);
	131	fOutput->Add(fHistSignal);
	132	fOutput->Add(fHistBackground);
	133
	134
8c34bd86	135	fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:Ptpi:PtK:Ptpi2:PtRec:PtTrue:PointingAngle:DecLeng:VxTrue:VxRec:VyRec,VzRec,InvMass:sigvert:d0Pi:d0K:d0Pi2");
	136	fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptkbkg:Ptpi2bkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:VyRecbkg:VzRecbkg:InvMassbkg:sigvertbkg:d0Pibkg:d0Kbkg:d0Pi2bkg");
	137
d486095a	138	return;
	139	}
	140
	141	//________________________________________________________________________
	142	void AliAnalysisTaskSEDplus::UserExec(Option_t /option*/)
	143	{
	144	// Execute analysis for current event:
	145	// heavy flavor candidates association to MC truth
	146
	147
	148
	149	AliAODEvent aod = dynamic_cast<AliAODEvent> (InputEvent());
	150
	151	// load Dplus->Kpipi candidates
	152	TClonesArray *array3Prong =
	153	(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
	154	if(!array3Prong) {
	155	printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
	156	return;
	157	}
	158
	159	// AOD primary vertex
	160	AliAODVertex vtx1 = (AliAODVertex)aod->GetPrimaryVertex();
	161	// vtx1->Print();
	162
	163	// load MC particles
	164	TClonesArray *arrayMC =
	165	(TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
	166	if(!arrayMC) {
	167	printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
	168	return;
	169	}
	170
	171	// load MC header
	172	AliAODMCHeader *mcHeader =
	173	(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
	174	if(!mcHeader) {
	175	printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
	176	return;
	177	}
	178	Int_t n3Prong = array3Prong->GetEntriesFast();
	179	printf("Number of D+->Kpipi: %d\n",n3Prong);
	180
	181
	182	for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
	183	AliAODRecoDecayHF3Prong d = (AliAODRecoDecayHF3Prong)array3Prong->UncheckedAt(i3Prong);
	184
	185
	186	Bool_t unsetvtx=kFALSE;
	187	if(!d->GetOwnPrimaryVtx()){
	188	d->SetOwnPrimaryVtx(vtx1);
	189	unsetvtx=kTRUE;
	190	}
8c34bd86	191	if(d->SelectDplus(fVHF->GetDplusCuts())) {
d486095a	192
d486095a	193
8c34bd86	194	Int_t labDp = d->MatchToMC(411,arrayMC);
	195
	196	if(labDp>=0) {
	197	AliAODMCParticle partDp = (AliAODMCParticle)arrayMC->At(labDp);
	198	Int_t pdgDp = TMath::Abs(partDp->GetPdgCode());
	199	if(pdgDp==411){
413e9336	200
8c34bd86	201	fHistSignal->Fill(d->InvMassDplus());
	202	fHistMass->Fill(d->InvMassDplus());
	203
	204	// Post the data already here
	205	PostData(1,fOutput);
	206
	207	AliAODMCParticle dg0 = (AliAODMCParticle)arrayMC->At(partDp->GetDaughter(0));
	208
	209	Float_t tmp[20];
	210	tmp[0]=pdgDp;
	211	tmp[1]=partDp->Px()-d->Px();
	212	tmp[2]=partDp->Py()-d->Py();
	213	tmp[3]=partDp->Pz()-d->Pz();
	214	tmp[4]=d->PtProng(0);
	215	tmp[5]=d->PtProng(1);
	216	tmp[6]=d->PtProng(2);
	217	tmp[7]=d->Pt();
	218	tmp[8]=partDp->Pt();
	219	tmp[9]=d->CosPointingAngle();
	220	tmp[10]=d->DecayLength();
	221	tmp[11]=dg0->Xv();
	222	tmp[12]=d->Xv();
	223	tmp[13]=d->Yv();
	224	tmp[14]=d->Zv();
	225	tmp[15]=d->InvMassDplus();
	226	tmp[16]=d->GetSigmaVert();
	227	tmp[17]=d->Getd0Prong(0);
	228	tmp[18]=d->Getd0Prong(1);
	229	tmp[19]=d->Getd0Prong(2);
	230
	231
	232	fNtupleDplus->Fill(tmp);
	233	PostData(2,fNtupleDplus);
	234	}
	235	} else {
	236
	237	Float_t tmpbkg[14];
	238	tmpbkg[0]=d->PtProng(0);
	239	tmpbkg[1]=d->PtProng(1);
	240	tmpbkg[2]=d->PtProng(2);
	241	tmpbkg[3]=d->Pt();
	242	tmpbkg[4]=d->CosPointingAngle();
	243	tmpbkg[5]=d->DecayLength();
	244	tmpbkg[6]=d->Xv();
	245	tmpbkg[7]=d->Yv();
	246	tmpbkg[8]=d->Zv();
	247	tmpbkg[9]=d->InvMassDplus();
	248	tmpbkg[10]=d->GetSigmaVert();
	249	tmpbkg[11]=d->Getd0Prong(0);
	250	tmpbkg[12]=d->Getd0Prong(1);
	251	tmpbkg[13]=d->Getd0Prong(2);
	252
	253	fHistBackground->Fill(d->InvMassDplus());
	254	fNtupleDplusbackg->Fill(tmpbkg);
	255	PostData(3,fNtupleDplusbackg);
413e9336	256	fHistMass->Fill(d->InvMassDplus());
fc8d975b	257	}
fc8d975b	258
8c34bd86	259
	260	}
	261
d486095a	262	if(unsetvtx) d->UnsetOwnPrimaryVtx();
	263
	264	}
	265	// end loop on D+->Kpipi
	266
	267
	268	return;
	269	}
	270
	271	//________________________________________________________________________
	272	void AliAnalysisTaskSEDplus::Terminate(Option_t /option*/)
	273	{
	274	// Terminate analysis
	275	//
	276	if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
	277
	278	fOutput = dynamic_cast<TList*> (GetOutputData(1));
	279	if (!fOutput) {
	280	printf("ERROR: fOutput not available\n");
	281	return;
	282	}
	283
d486095a	284	fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
	285	fHistSignal = dynamic_cast<TH1F*>(fOutput->FindObject("fHistSignal"));
	286	fHistBackground = dynamic_cast<TH1F*>(fOutput->FindObject("fHistBackground"));
8c34bd86	287	fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
8c34bd86	288	fNtupleDplusbackg = dynamic_cast<TNtuple*>(GetOutputData(3));
fc8d975b	289
8c34bd86	290	return;
d486095a	291	}
d486095a	292