[u/mrichter/AliRoot.git] / HLT / trigger / AliD0Trigger.C

/****************************************************************************
 *                                                                          *
 * This macro computes the position of the D0 decay vertex using            *
 * helix DCA minimization by J.Belikov                                      *
 *                                                                          *
 * Reconstructed D0 are stored in a tree as AliD0toKpi objects              *
 *                                                                          *
 ****************************************************************************/ 
//#define __COMPILE__
//#ifdef __COMPILE__
#if !defined(__CINT__) || defined(__MAKECINT__)
//-- --- standard headers------------- 
#include <iostream.h>
//--------Root headers ---------------
#include <TSystem.h>
#include <TFile.h>
#include <TString.h>
#include <TStopwatch.h>
#include <TObject.h>
#include <TVector3.h>
#include <TTree.h>
#include <TParticle.h>
#include <TArray.h>
//----- AliRoot headers ---------------
#include "alles.h"
#include "AliRun.h"
#include "AliKalmanTrack.h"
#include "AliITStrackV2.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliV0vertex.h"
#include "AliV0vertexer.h"
#include "AliITSVertex.h"
#include "AliITSVertexer.h"
#include "AliITSVertexerTracks.h"
#include "AliD0toKpi.h"
#endif
//-------------------------------------

typedef struct {
  Int_t lab;
  Int_t pdg;
  Int_t mumlab;
  Int_t mumpdg;
  Float_t Vx,Vy,Vz;
  Float_t Px,Py,Pz;
} RECTRACK;

// field (T)
const Double_t kBz = 0.4;

// primary vertex
Double_t gv1[3] = {0.,0.,0,};

// sigle track cuts
const Double_t kPtCut = 0.5;  // GeV/c
const Double_t kd0Cut = 50.; // micron
//const Double_t kPtCut = 0.;  // GeV/c
//const Double_t kd0Cut = 0.; // micron

// cuts on D0 candidate (to be passed to function AliD0toKpi::Select())
//
// cuts[0] = inv. mass half width [GeV]   
// cuts[1] = dca [micron]
// cuts[2] = cosThetaStar 
// cuts[3] = pTK [GeV/c]
// cuts[4] = pTPi [GeV/c]
// cuts[5] = d0K [micron]   upper limit!
// cuts[6] = d0Pi [micron]  upper limit!
// cuts[7] = d0d0 [micron^2]
// cuts[8] = cosThetaPoint
const Double_t kD0cuts[9] = {0.012,
                             500.,
			     0.8,
			     0.5,
			     0.5,
			     500.,
			     500.,
			    -5000.,
			     0.8};
const Double_t kD0cutsAll[9] = {.1,
				1.e5,
				1.1,
				0.,
				0.,
				1.e10,
				1.e10,
				1.e10,
				-1.1};
// mass cut for 1st level rejection
const Double_t kMassCut = 0.1; // GeV

// this function ckecks files existence
Bool_t GetInputFiles(char* path[1024]);

// 1st level rejection based on inv. mass
Bool_t SelectInvMass(const Double_t p[6]);

// this function applies single track cuts
Bool_t TrkCuts(const AliITStrackV2& trk);

// this function creates TObjArrays with positive and negative tracks
void   SelectTracks(TTree& itsTree,
		      TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
		      TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);


void   AliD0Trigger(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {

  const Char_t *name="AliD0Trigger";
  cerr<<'\n'<<name<<"...\n";
  gBenchmark->Start(name);  
  TStopwatch timer;
  AliKalmanTrack::SetConvConst(100/0.299792458/kBz);

  // load library for D0 candidates class
  //gSystem->Load("libD0toKpi.so");

  // check existence of input files 
  if(!GetInputFiles(path)) { cerr<<"No tracks file found"<<endl; return; }

  Char_t fname[1024];
  sprintf(fname,"%s/AliD0toKpiBkg.root",path);
  TString outName(fname);

  Int_t ev;
  Int_t nD0rec=0,nD0rec1ev=0;
  Double_t dca;
  Double_t v2[3],mom[6],d0[2];
  Int_t pdg[2],mum[2];
  Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
  Int_t iTrkP,iTrkN,itsEntries;
  Int_t nTrksP=0,nTrksN=0;
  Int_t okD0=0,okD0bar=0;
  Char_t trksName[100],refsName[100],vtxName[100];
  RECTRACK rectrk;  
  AliITStrackV2 *postrack = 0;
  AliITStrackV2 *negtrack = 0;

  // create the AliITSVertexerTracks object
  AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
  //vertexer1->SetUseThrustFrame(1);
  vertexer1->SetMinTracks(2);
  vertexer1->SetDebug(0);
  Int_t skipped[2];

  // define the cuts for vertexing
  Double_t vtxcuts[]={33., // max. allowed chi2
		      0.0, // min. allowed negative daughter's impact param 
		      0.0, // min. allowed positive daughter's impact param 
		      1.0, // max. allowed DCA between the daughter tracks
        	     -1.0, // min. allowed cosine of V0's pointing angle
		      0.0, // min. radius of the fiducial volume
		      2.9};// max. radius of the fiducial volume
  
  // create the AliV0vertexer object
  AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);


  // Create tree for reconstructed D0s
  TTree D0Tree("TreeD0","Tree with D0 candidates");
  AliD0toKpi *ioD0toKpi=0;
  D0Tree.Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);

  // Open file with ITS tracks
  Char_t fnameTrack[1024];
  sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
  TFile* itstrks = TFile::Open(fnameTrack);

  // Open file with ITS track references
  Char_t fnameRef[1024];
  sprintf(fnameRef,"%s/ITStracksRefFile.root",path);
  TFile* itsrefs = TFile::Open(fnameRef);

  // Open Primary Vertex File
  Char_t fnameVertex[1024];
  sprintf(fnameVertex,"%s/Vtx4Tracking.root",path);
  TFile* vertexFile = TFile::Open(fnameVertex);

  // loop on events in file
  for(ev=evFirst; ev<=evLast; ev++) {
    printf(" --- Processing event  %d ---\n",ev);
    sprintf(vtxName,"VertexTracks_%d",ev);
    sprintf(trksName,"TreeT_ITS_%d",ev);
    sprintf(refsName,"Tree_Ref_%d",ev);


    // tracks from ITS
    TTree *itsTree=(TTree*)itstrks->Get(trksName);
    if(!itsTree) continue;
    itsEntries = (Int_t)itsTree->GetEntries();
    printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
    
    // tree from reference file
    TTree *refTree=(TTree*)itsrefs->Get(refsName);
    refTree->SetBranchAddress("rectracks",&rectrk);
    
    //get the primary vertex from headerfile

    AliITSVertex *vertex = vertexFile->Get("Vertex_0");
    vertex->GetXYZ(gv1);

    cout<<gv1[2]<<endl;
    
    // call function which applies sigle track selection and
    // separetes positives and negatives
    
    TObjArray trksP(itsEntries/2);
    Int_t *itsEntryP = new Int_t[itsEntries];
    TObjArray trksN(itsEntries/2);
    Int_t *itsEntryN = new Int_t[itsEntries];
    SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);

    nD0rec1ev = 0;

    cout<<"NnegTracks: "<<nTrksN<<endl;

    // loop on positive tracks
    for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
      //if(iTrkP % 10==0) cerr<<"--- Processing positive track number: "<<iTrkP<<" of "<<nTrksP<<"\r";
	  
      // get track from track array
      postrack = (AliITStrackV2*)trksP.At(iTrkP);

      // get info on tracks PDG and mothers PDG from reference file
      refTree->GetEvent(itsEntryP[iTrkP]);
      pdg[0] = rectrk.pdg;
      mum[0] = rectrk.mumpdg;

      // loop on negative tracks 
      for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
	
	// get track from track array
	negtrack = (AliITStrackV2*)trksN.At(iTrkN);
      
	// get info on tracks PDG and mothers PDG from reference file
	refTree->GetEvent(itsEntryN[iTrkN]);
	pdg[1] = rectrk.pdg;
	mum[1] = rectrk.mumpdg;
 
	//
	// ----------- DCA MINIMIZATION ------------------
	//
	// find the DCA and propagate the tracks to the DCA 
	dca = vertexer2->PropagateToDCA(negtrack,postrack);

	// define the AliV0vertex object
	AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
	  
	// get position of the vertex
	vertex2->GetXYZ(v2[0],v2[1],v2[2]);
	
	delete vertex2;
  
	// momenta of the tracks at the vertex
	ptP = 1./TMath::Abs(postrack->Get1Pt());
	alphaP = postrack->GetAlpha();
	phiP = alphaP+TMath::ASin(postrack->GetSnp());
	mom[0] = ptP*TMath::Cos(phiP); 
	mom[1] = ptP*TMath::Sin(phiP);
	mom[2] = ptP*postrack->GetTgl();
	  
	ptN = 1./TMath::Abs(negtrack->Get1Pt());
	alphaN = negtrack->GetAlpha();
	phiN = alphaN+TMath::ASin(negtrack->GetSnp());
	mom[3] = ptN*TMath::Cos(phiN); 
	mom[4] = ptN*TMath::Sin(phiN);
	mom[5] = ptN*negtrack->GetTgl();
	  

	Bool_t goodVtx = kFALSE;
	// no vertexing if DeltaMass > kMassCut 
	if(SelectInvMass(mom)) {
	  // primary vertex from other tracks in event
	  vertexer1->SetVtxStart(0.,0.);
	  skipped[0] = itsEntryP[iTrkP];
	  skipped[1] = itsEntryN[iTrkN];
	  vertexer1->SetSkipTracks(2,skipped);
	  // impact parameters of the tracks w.r.t. the smeared primary vertex
	  d0[0] =  10000.*postrack->GetD(gv1[0],gv1[1]);
	  d0[1] = -10000.*negtrack->GetD(gv1[0],gv1[1]);
	}
	  
	// create the object TD0rec and store it in the tree
	AliD0toKpi theD0(kFALSE,ev,gv1,v2,dca,mom,d0,pdg,mum);

	
	// select D0s
	if(goodVtx && theD0.Select(kD0cuts,okD0,okD0bar)) {
	      
	  // compute the weights
	  theD0.ComputeWgts();

	  // fill the tree
	  ioD0toKpi=&theD0;
	  D0Tree.Fill();

	  nD0rec++; nD0rec1ev++;
	  ioD0toKpi=0;  
	}
	
	negtrack = 0;
      } // loop on negative tracks
      postrack = 0;
    }   // loop on positive tracks

    delete [] itsEntryP;
    delete [] itsEntryN;
    delete itsTree;
    delete refTree;

    printf("\n+++\n+++ Number of D0 candidates: %d\n+++\n",nD0rec1ev);
    printf("----------------------------------------------------------\n");

  }    // loop on events in file

  printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);

  delete vertexer1;
  delete vertexer2;

  itstrks->Close();
  itsrefs->Close();

  // store tree in a file
  TFile* outroot = new TFile(outName.Data(),"recreate");
  D0Tree.Write();
  outroot->Close();
  delete outroot;

  gBenchmark->Stop(name);
  gBenchmark->Show(name);

  return;
}
//___________________________________________________________________________
Bool_t   GetInputFiles(Char_t* path2="./") {
  Char_t fnameTrack[1024];
  sprintf(fnameTrack,"%s/AliITStracksV2.root",path2);
  TString itsName(fnameTrack);
  Char_t fnameRef[1024];
  sprintf(fnameRef,"%s/ITStracksRefFile.root",path2);
  TString refName(fnameRef);

  if(gSystem->AccessPathName(itsName.Data(),kFileExists) ||
     gSystem->AccessPathName(refName.Data(),kFileExists)) return kFALSE;

  return kTRUE;
}
//___________________________________________________________________________
Bool_t SelectInvMass(const Double_t p[6]) {

  Double_t mD0 = 1.8645;
  Double_t mPi = 0.13957;
  Double_t mKa = 0.49368;

  Double_t energy[2];
  Double_t mom2[2],momTot2;

  mom2[0] = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
  mom2[1] = p[3]*p[3]+p[4]*p[4]+p[5]*p[5];

  momTot2 = (p[0]+p[3])*(p[0]+p[3])+
            (p[1]+p[4])*(p[1]+p[4])+
            (p[2]+p[5])*(p[2]+p[5]);

  // D0 -> K- Pi+
  energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
  energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);

  Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
    

  // D0bar -> K+ Pi-
  energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
  energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);

  Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);

  if(TMath::Abs(minvD0-mD0) < kMassCut)    return kTRUE;
  if(TMath::Abs(minvD0bar-mD0) < kMassCut) return kTRUE;
  return kFALSE;
}
//___________________________________________________________________________
void   SelectTracks(TTree& itsTree,
		    TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
		    TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
//
// this function creates two TObjArrays with positive and negative tracks
//
  nTrksP=0,nTrksN=0;

 
  Int_t entr = (Int_t)itsTree.GetEntries();

  // trasfer tracks from tree to arrays
  for(Int_t i=0; i<entr; i++) {

    AliITStrackV2 *itstrack = new AliITStrackV2; 
    itsTree.SetBranchAddress("tracks",&itstrack);

    itsTree.GetEvent(i);

    // single track selection
    if(!TrkCuts(*itstrack)) { delete itstrack; continue; }

    if(itstrack->Get1Pt()>0.) { // negative track
      trksN.AddLast(itstrack);
      itsEntryN[nTrksN] = i;
      nTrksN++;
    } else {                    // positive track
      trksP.AddLast(itstrack);
      itsEntryP[nTrksP] = i;
      nTrksP++;
    }

  }

  return;
}
//____________________________________________________________________________
Bool_t TrkCuts(const AliITStrackV2& trk) {
// 
// this function tells if track passes some kinematical cuts  
//
  if(TMath::Abs(1./trk.Get1Pt()) < kPtCut)                return kFALSE;
  if(TMath::Abs(10000.*trk.GetD(gv1[0],gv1[1])) < kd0Cut) return kFALSE;

  return kTRUE;
}
Commit	Line	Data
6f388e0d	1	/****************************************************************************
	2	* *
	3	* This macro computes the position of the D0 decay vertex using *
	4	* helix DCA minimization by J.Belikov *
	5	* *
	6	* Reconstructed D0 are stored in a tree as AliD0toKpi objects *
	7	* *
	8	****************************************************************************/
	9	//#define __COMPILE__
	10	//#ifdef __COMPILE__
	11	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	12	//-- --- standard headers-------------
	13	#include <iostream.h>
	14	//--------Root headers ---------------
	15	#include <TSystem.h>
	16	#include <TFile.h>
	17	#include <TString.h>
	18	#include <TStopwatch.h>
	19	#include <TObject.h>
	20	#include <TVector3.h>
	21	#include <TTree.h>
	22	#include <TParticle.h>
	23	#include <TArray.h>
	24	//----- AliRoot headers ---------------
	25	#include "alles.h"
	26	#include "AliRun.h"
	27	#include "AliKalmanTrack.h"
	28	#include "AliITStrackV2.h"
	29	#include "AliHeader.h"
	30	#include "AliGenEventHeader.h"
	31	#include "AliV0vertex.h"
	32	#include "AliV0vertexer.h"
	33	#include "AliITSVertex.h"
	34	#include "AliITSVertexer.h"
	35	#include "AliITSVertexerTracks.h"
	36	#include "AliD0toKpi.h"
	37	#endif
	38	//-------------------------------------
	39
	40	typedef struct {
	41	Int_t lab;
	42	Int_t pdg;
	43	Int_t mumlab;
	44	Int_t mumpdg;
	45	Float_t Vx,Vy,Vz;
	46	Float_t Px,Py,Pz;
	47	} RECTRACK;
	48
	49	// field (T)
	50	const Double_t kBz = 0.4;
	51
	52	// primary vertex
	53	Double_t gv1[3] = {0.,0.,0,};
	54
	55	// sigle track cuts
	56	const Double_t kPtCut = 0.5; // GeV/c
	57	const Double_t kd0Cut = 50.; // micron
	58	//const Double_t kPtCut = 0.; // GeV/c
	59	//const Double_t kd0Cut = 0.; // micron
	60
	61	// cuts on D0 candidate (to be passed to function AliD0toKpi::Select())
	62	//
	63	// cuts[0] = inv. mass half width [GeV]
	64	// cuts[1] = dca [micron]
65	// cuts[2] = cosThetaStar
66	// cuts[3] = pTK [GeV/c]
67	// cuts[4] = pTPi [GeV/c]
68	// cuts[5] = d0K [micron] upper limit!
69	// cuts[6] = d0Pi [micron] upper limit!
70	// cuts[7] = d0d0 [micron^2]
71	// cuts[8] = cosThetaPoint
72	const Double_t kD0cuts[9] = {0.012,
73	500.,
74	0.8,
75	0.5,
76	0.5,
77	500.,
78	500.,
79	-5000.,
80	0.8};
81	const Double_t kD0cutsAll[9] = {.1,
82	1.e5,
83	1.1,
84	0.,
85	0.,
86	1.e10,
87	1.e10,
88	1.e10,
89	-1.1};
90	// mass cut for 1st level rejection
91	const Double_t kMassCut = 0.1; // GeV
92
93	// this function ckecks files existence
94	Bool_t GetInputFiles(char* path[1024]);
95
96	// 1st level rejection based on inv. mass
97	Bool_t SelectInvMass(const Double_t p[6]);
98
99	// this function applies single track cuts
100	Bool_t TrkCuts(const AliITStrackV2& trk);
101
102	// this function creates TObjArrays with positive and negative tracks
103	void SelectTracks(TTree& itsTree,
104	TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
105	TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
106
107
108	void AliD0Trigger(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {
109
110	const Char_t *name="AliD0Trigger";
111	cerr<<'\n'<<name<<"...\n";
112	gBenchmark->Start(name);
113	TStopwatch timer;
114	AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
115
116	// load library for D0 candidates class
117	//gSystem->Load("libD0toKpi.so");
118
119	// check existence of input files
120	if(!GetInputFiles(path)) { cerr<<"No tracks file found"<<endl; return; }
121
122	Char_t fname[1024];
123	sprintf(fname,"%s/AliD0toKpiBkg.root",path);
124	TString outName(fname);
125
126	Int_t ev;
127	Int_t nD0rec=0,nD0rec1ev=0;
128	Double_t dca;
129	Double_t v2[3],mom[6],d0[2];
130	Int_t pdg[2],mum[2];
131	Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
132	Int_t iTrkP,iTrkN,itsEntries;
133	Int_t nTrksP=0,nTrksN=0;
134	Int_t okD0=0,okD0bar=0;
135	Char_t trksName[100],refsName[100],vtxName[100];
136	RECTRACK rectrk;
137	AliITStrackV2 *postrack = 0;
138	AliITStrackV2 *negtrack = 0;
139
140	// create the AliITSVertexerTracks object
141	AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
142	//vertexer1->SetUseThrustFrame(1);
143	vertexer1->SetMinTracks(2);
144	vertexer1->SetDebug(0);
145	Int_t skipped[2];
146
147	// define the cuts for vertexing
148	Double_t vtxcuts[]={33., // max. allowed chi2
149	0.0, // min. allowed negative daughter's impact param
150	0.0, // min. allowed positive daughter's impact param
151	1.0, // max. allowed DCA between the daughter tracks
152	-1.0, // min. allowed cosine of V0's pointing angle
153	0.0, // min. radius of the fiducial volume
154	2.9};// max. radius of the fiducial volume
155
156	// create the AliV0vertexer object
157	AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
158
159
160	// Create tree for reconstructed D0s
161	TTree D0Tree("TreeD0","Tree with D0 candidates");
162	AliD0toKpi *ioD0toKpi=0;
163	D0Tree.Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
164
165	// Open file with ITS tracks
166	Char_t fnameTrack[1024];
167	sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
168	TFile* itstrks = TFile::Open(fnameTrack);
169
170	// Open file with ITS track references
171	Char_t fnameRef[1024];
172	sprintf(fnameRef,"%s/ITStracksRefFile.root",path);
173	TFile* itsrefs = TFile::Open(fnameRef);
174
175	// Open Primary Vertex File
176	Char_t fnameVertex[1024];
177	sprintf(fnameVertex,"%s/Vtx4Tracking.root",path);
178	TFile* vertexFile = TFile::Open(fnameVertex);
179
180	// loop on events in file
181	for(ev=evFirst; ev<=evLast; ev++) {
182	printf(" --- Processing event %d ---\n",ev);
183	sprintf(vtxName,"VertexTracks_%d",ev);
184	sprintf(trksName,"TreeT_ITS_%d",ev);
185	sprintf(refsName,"Tree_Ref_%d",ev);
186
187
188	// tracks from ITS
189	TTree itsTree=(TTree)itstrks->Get(trksName);
190	if(!itsTree) continue;
191	itsEntries = (Int_t)itsTree->GetEntries();
192	printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
193
194	// tree from reference file
195	TTree refTree=(TTree)itsrefs->Get(refsName);
196	refTree->SetBranchAddress("rectracks",&rectrk);
197
198	//get the primary vertex from headerfile
199
200	AliITSVertex *vertex = vertexFile->Get("Vertex_0");
201	vertex->GetXYZ(gv1);
202
203	cout<<gv1[2]<<endl;
204
205	// call function which applies sigle track selection and
206	// separetes positives and negatives
207
208	TObjArray trksP(itsEntries/2);
209	Int_t *itsEntryP = new Int_t[itsEntries];
210	TObjArray trksN(itsEntries/2);
211	Int_t *itsEntryN = new Int_t[itsEntries];
212	SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
213
214	nD0rec1ev = 0;
215
216	cout<<"NnegTracks: "<<nTrksN<<endl;
217
218	// loop on positive tracks
219	for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
220	//if(iTrkP % 10==0) cerr<<"--- Processing positive track number: "<<iTrkP<<" of "<<nTrksP<<"\r";
221
222	// get track from track array
223	postrack = (AliITStrackV2*)trksP.At(iTrkP);
224
225	// get info on tracks PDG and mothers PDG from reference file
226	refTree->GetEvent(itsEntryP[iTrkP]);
227	pdg[0] = rectrk.pdg;
228	mum[0] = rectrk.mumpdg;
229
230	// loop on negative tracks
231	for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
232
233	// get track from track array
234	negtrack = (AliITStrackV2*)trksN.At(iTrkN);
235
236	// get info on tracks PDG and mothers PDG from reference file
237	refTree->GetEvent(itsEntryN[iTrkN]);
238	pdg[1] = rectrk.pdg;
239	mum[1] = rectrk.mumpdg;
240
241	//
242	// ----------- DCA MINIMIZATION ------------------
243	//
244	// find the DCA and propagate the tracks to the DCA
245	dca = vertexer2->PropagateToDCA(negtrack,postrack);
246
247	// define the AliV0vertex object
248	AliV0vertex vertex2 = new AliV0vertex(negtrack,*postrack);
249
250	// get position of the vertex
251	vertex2->GetXYZ(v2[0],v2[1],v2[2]);
252
253	delete vertex2;
254
255	// momenta of the tracks at the vertex
256	ptP = 1./TMath::Abs(postrack->Get1Pt());
257	alphaP = postrack->GetAlpha();
258	phiP = alphaP+TMath::ASin(postrack->GetSnp());
259	mom[0] = ptP*TMath::Cos(phiP);
260	mom[1] = ptP*TMath::Sin(phiP);
261	mom[2] = ptP*postrack->GetTgl();
262
263	ptN = 1./TMath::Abs(negtrack->Get1Pt());
264	alphaN = negtrack->GetAlpha();
265	phiN = alphaN+TMath::ASin(negtrack->GetSnp());
266	mom[3] = ptN*TMath::Cos(phiN);
267	mom[4] = ptN*TMath::Sin(phiN);
268	mom[5] = ptN*negtrack->GetTgl();
269
270
271	Bool_t goodVtx = kFALSE;
272	// no vertexing if DeltaMass > kMassCut
273	if(SelectInvMass(mom)) {
274	// primary vertex from other tracks in event
275	vertexer1->SetVtxStart(0.,0.);
276	skipped[0] = itsEntryP[iTrkP];
277	skipped[1] = itsEntryN[iTrkN];
278	vertexer1->SetSkipTracks(2,skipped);
279	// impact parameters of the tracks w.r.t. the smeared primary vertex
280	d0[0] = 10000.*postrack->GetD(gv1[0],gv1[1]);
281	d0[1] = -10000.*negtrack->GetD(gv1[0],gv1[1]);
282	}
283
284	// create the object TD0rec and store it in the tree
285	AliD0toKpi theD0(kFALSE,ev,gv1,v2,dca,mom,d0,pdg,mum);
286
287
288	// select D0s
289	if(goodVtx && theD0.Select(kD0cuts,okD0,okD0bar)) {
290
291	// compute the weights
292	theD0.ComputeWgts();
293
294	// fill the tree
295	ioD0toKpi=&theD0;
296	D0Tree.Fill();
297
298	nD0rec++; nD0rec1ev++;
299	ioD0toKpi=0;
300	}
301
302	negtrack = 0;
303	} // loop on negative tracks
304	postrack = 0;
305	} // loop on positive tracks
306
307	delete [] itsEntryP;
308	delete [] itsEntryN;
309	delete itsTree;
310	delete refTree;
311
312	printf("\n+++\n+++ Number of D0 candidates: %d\n+++\n",nD0rec1ev);
313	printf("----------------------------------------------------------\n");
314
315	} // loop on events in file
316
317	printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
318
319	delete vertexer1;
320	delete vertexer2;
321
322	itstrks->Close();
323	itsrefs->Close();
324
325	// store tree in a file
326	TFile* outroot = new TFile(outName.Data(),"recreate");
327	D0Tree.Write();
328	outroot->Close();
329	delete outroot;
330
331	gBenchmark->Stop(name);
332	gBenchmark->Show(name);
333
334	return;
335	}
336	//___________________________________________________________________________
337	Bool_t GetInputFiles(Char_t* path2="./") {
338	Char_t fnameTrack[1024];
339	sprintf(fnameTrack,"%s/AliITStracksV2.root",path2);
340	TString itsName(fnameTrack);
341	Char_t fnameRef[1024];
342	sprintf(fnameRef,"%s/ITStracksRefFile.root",path2);
343	TString refName(fnameRef);
344
345	if(gSystem->AccessPathName(itsName.Data(),kFileExists) \|\|
346	gSystem->AccessPathName(refName.Data(),kFileExists)) return kFALSE;
347
348	return kTRUE;
349	}
350	//___________________________________________________________________________
351	Bool_t SelectInvMass(const Double_t p[6]) {
352
353	Double_t mD0 = 1.8645;
354	Double_t mPi = 0.13957;
355	Double_t mKa = 0.49368;
356
357	Double_t energy[2];
358	Double_t mom2[2],momTot2;
359
360	mom2[0] = p[0]p[0]+p[1]p[1]+p[2]*p[2];
361	mom2[1] = p[3]p[3]+p[4]p[4]+p[5]*p[5];
362
363	momTot2 = (p[0]+p[3])*(p[0]+p[3])+
364	(p[1]+p[4])*(p[1]+p[4])+
365	(p[2]+p[5])*(p[2]+p[5]);
366
367	// D0 -> K- Pi+
368	energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
369	energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
370
371	Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
372
373
374	// D0bar -> K+ Pi-
375	energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
376	energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
377
378	Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
379
380	if(TMath::Abs(minvD0-mD0) < kMassCut) return kTRUE;
381	if(TMath::Abs(minvD0bar-mD0) < kMassCut) return kTRUE;
382	return kFALSE;
383	}
384	//___________________________________________________________________________
385	void SelectTracks(TTree& itsTree,
386	TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
387	TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
388	//
389	// this function creates two TObjArrays with positive and negative tracks
390	//
391	nTrksP=0,nTrksN=0;
392
393
394	Int_t entr = (Int_t)itsTree.GetEntries();
395
396	// trasfer tracks from tree to arrays
397	for(Int_t i=0; i<entr; i++) {
398
399	AliITStrackV2 *itstrack = new AliITStrackV2;
400	itsTree.SetBranchAddress("tracks",&itstrack);
401
402	itsTree.GetEvent(i);
403
404	// single track selection
405	if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
406
407	if(itstrack->Get1Pt()>0.) { // negative track
408	trksN.AddLast(itstrack);
409	itsEntryN[nTrksN] = i;
410	nTrksN++;
411	} else { // positive track
412	trksP.AddLast(itstrack);
413	itsEntryP[nTrksP] = i;
414	nTrksP++;
415	}
416
417	}
418
419	return;
420	}
421	//____________________________________________________________________________
422	Bool_t TrkCuts(const AliITStrackV2& trk) {
423	//
424	// this function tells if track passes some kinematical cuts
425	//
426	if(TMath::Abs(1./trk.Get1Pt()) < kPtCut) return kFALSE;
427	if(TMath::Abs(10000.*trk.GetD(gv1[0],gv1[1])) < kd0Cut) return kFALSE;
428
429	return kTRUE;
430	}
431
432
433
434
435
436
437
438