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

//#define __COMPILE__
//#ifdef __COMPILE__
#if !defined(__CINT__) || defined(__MAKECINT__)
//-- --- standard headers------------- 
#include <iostream.h>
#include <fstream.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 <TObjArray>
#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 "AliD0Trigger.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 primaryvertex[3]={0.,0.,0.};

//sec. vertex
double v2[3]={0,0,0};

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


//cuts for combined tracks
const Double_t cuts[7] = {0.005,     // 0.005 cuts[0] = lowest V0 cut  (cm)
			  0.800,     // 0.015 cuts[1] = highest V0 cut (cm)
			  0.012,     // 0.012 cuts[2] = inv. mass cut (diferense) (Gev/c)
			  0.8,       // 0.8   cuts[3] = min. cosine value for pointing angle
			  -5000,     // -5000 cuts[4] = d0d0
			  0,         // 0.8   cuts[5] = costhetastar
			  0.5};      // 0.5   cuts[6] = ptchild
//cut for distance of closest aprach
double cutDCA=0.05;   //0.05

// 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 GetPrimaryVertex(int i,Char_t* path="./");

//void PtD0(Char_t* path="./");

Int_t iTrkP,iTrkN,itsEntries;
Char_t trksName[100],refsName[100];
Int_t nTrksP=0,nTrksN=0;
Int_t nD0=0;
int ev=0;
double mom[6];
int event[10000];
int index=0;
Bool_t isSignal;
Int_t nTotEv=0,nD0recSgn=0,nD0recBkgS=0,nD0recBkg=0,nD0rec1ev=0;
Int_t pdg[2],mum[2],mumlab[2];
RECTRACK rectrk;

void RunD0offline(Char_t* path="./",bool h=false,bool PtD0=false) {
    
  AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
  
  // Open file with ITS tracks
  Char_t fnameTrack[1024];
  //sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
  sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
  TFile* itstrks = TFile::Open(fnameTrack);
  
  Char_t refFile[1024];
  //sprintf(fnameTrack,"%s/ITStracksRefFile.root",path);
  sprintf(refFile,"%s/ITStracksRefFile.root",path);
  TFile* itsrefs = TFile::Open(refFile);
  

  //the offline stuff
  // 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
  
  TH1F *h1 = new TH1F("h1","Transvers momentun of reconstructed D0",100,0,10);
  TH1F *h2 = new TH1F("h2","Transvers momentun of D0 with |Eta|<0.9",100,0,10);
  TH1F *h3 = new TH1F("h3","Eta reconstructed of D0",100,-5,5);
  TH1F *h4 = new TH1F("h4","Eta of D0",100,-5,5);
  
  Char_t falice[1024];
  sprintf(falice,"%s/galice.root",path);
  TFile *f = new TFile(falice);   
  gAlice=(AliRun*)f->Get("gAlice");
  int nEvent=gAlice->GetEventsPerRun();
  //int nEvent=20;
  cout<<"#Events: "<<nEvent<<endl;
  delete gAlice;

  const Char_t *name="AliD0offline";
  cerr<<'\n'<<name<<"...\n";
  gBenchmark->Start(name); 
  
  for(ev=0;ev<nEvent;ev++) {
    
    cout<<"\n Event: "<<ev<<endl;
    
    // tracks from ITS
    sprintf(trksName,"TreeT_ITS_%d",ev);
    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
    
    sprintf(refsName,"Tree_Ref_%d",ev);
    TTree *refTree=(TTree*)itsrefs->Get(refsName);
    refTree->SetBranchAddress("rectracks",&rectrk);

    //Getting primary Vertex
    GetPrimaryVertex(ev,path);
    
    // count the total number of events
    nTotEv++;

    // 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); 
    
    cout<<"#pos: "<<nTrksP<<endl;
    cout<<"#neg: "<<nTrksN<<endl;
    
    // create the AliV0vertexer object
    AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
    
    AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
    
    double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
    
    for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
      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;
      mumlab[0] = rectrk.mumlab;

      for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
	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;
	mumlab[1] = rectrk.mumlab;

	D0->SetTracks(postrack,negtrack);

	// ----------- DCA MINIMIZATION ------------------
	//
	// find the DCA and propagate the tracks to the DCA 
	double dca = vertexer2->PropagateToDCA(negtrack,postrack);
	
	if(dca<cutDCA){
	  // 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;	
	  //if(D0->FindV0offline(v2) && D0->d0d0()){
	  if(D0->d0d0()){
	    D0->SetV0(v2);
	    D0->FindMomentaOffline();
	    if(D0->FindInvMass() && D0->CosThetaStar() && D0->PointingAngle() && D0->pTchild()){
	      nD0++;
	      event[index]=ev; index++;
	      if(h){
		h1->Fill(D0->Pt());
		h3->Fill(D0->Eta());
	      }
	      
	      if(mumlab[0]==mumlab[1] && TMath::Abs(mum[0])==421 && TMath::Abs(mum[1])==421) { 
		nD0recSgn++;
	      } 
	      else { 
		if(TMath::Abs(mum[0])==421 || TMath::Abs(mum[0])==411 ||
		   TMath::Abs(mum[1])==421 || TMath::Abs(mum[1])==411) {
		  nD0recBkgS++;
		} else {
		  nD0recBkg++;
		}  
	      }
	    }
	  } 
	}
      }
    } 
    
    //delete D0;
    //delete itstrks;
    //delete itsTree;
    //delete trksP;
    //delete itsEntryP;
    //delete trksN;
    //delete itsEntryN;
  }
  
  cout<<"\nMy #D0: "<<nD0<<"\n"<<endl;
  printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
  printf("\n+++\n+++ Total number of D0 candidates:\n
                 +++    Sgn:   %d\n
                 +++    BkgS:  %d\n
                 +++    Bkg:   %d\n+++\n",nD0recSgn,nD0recBkgS,nD0recBkg);

  gBenchmark->Stop(name);
  gBenchmark->Show(name);
  
  if(h){
    if(!PtD0){
      TCanvas *c = new TCanvas("c","c",700,1000);
      c->Divide(1,2);
      c->cd(1);
      h1->Draw();
      pt = new TPaveText(7.3,0.4,11,1.8,"br");
      pt->SetFillColor(0);
      pt->SetBorderSize(1);
      pt->AddText("Cuts:");
      Char_t st[1024];
      sprintf(st,"First Pt:         %g",kPtCut);
      pt->AddText(st);
      sprintf(st,"d0 low:           %g",kd0Cut);
      pt->AddText(st);
      sprintf(st,"d0 high:          %g",kd0CutHigh);
      pt->AddText(st);
      sprintf(st,"V0 low:           %g",cuts[0]);
      pt->AddText(st);
      sprintf(st,"V0 high:          %g",cuts[1]);
      pt->AddText(st);
      sprintf(st,"InvMass Diff:     %g",cuts[2]);
      pt->AddText(st);
      sprintf(st,"cosPointingAngle: %g",cuts[3]);
      pt->AddText(st);
      sprintf(st,"d0d0:             %g",cuts[4]);
      pt->AddText(st);
      sprintf(st,"cosThetaStar:     %g",cuts[5]);
      pt->AddText(st);
      sprintf(st,"PtChild:          %g",cuts[6]);
      pt->AddText(st);
      sprintf(st,"DCA:              %g",cutDCA);
      pt->AddText(st);
      pt->Draw();
      c_1->Modified();
      c->cd();
      c->cd(2);
      h3->Draw();
    }
    if(PtD0){
      PtD0(path,h2,h4);
      TCanvas *c = new TCanvas("c","c",1000,700);
      c->Divide(2,2);
      c->cd(1);
      h1->Draw();
      pt = new TPaveText(7.3,0.4,11,1.8,"br");
      pt->SetFillColor(0);
      pt->SetBorderSize(1);
      pt->AddText("Cuts:");
      Char_t st[1024];
      sprintf(st,"First Pt:         %g",kPtCut);
      pt->AddText(st);
      sprintf(st,"d0 low:           %g",kd0Cut);
      pt->AddText(st);
      sprintf(st,"d0 high:          %g",kd0CutHigh);
      pt->AddText(st);
      sprintf(st,"V0 low:           %g",cuts[0]);
      pt->AddText(st);
      sprintf(st,"V0 high:          %g",cuts[1]);
      pt->AddText(st);
      sprintf(st,"InvMass Diff:     %g",cuts[2]);
      pt->AddText(st);
      sprintf(st,"cosPointAng: %g",cuts[3]);
      pt->AddText(st);
      sprintf(st,"d0d0:             %g",cuts[4]);
      pt->AddText(st);
      sprintf(st,"cosTheta*:     %g",cuts[5]);
      pt->AddText(st);
      sprintf(st,"PtChild:          %g",cuts[6]);
      pt->AddText(st);
      sprintf(st,"DCA:              %g",cutDCA);
      pt->AddText(st);
      pt->Draw();
      c_1->Modified();
      c->cd();
      c->cd(2);
      h3->Draw();
      c->cd(3);
      h2->Draw();
      c->cd(4);
      h4->Draw();
    }    
  }
  if(h){
    if(!PtD0){
      Char_t outName[1024];
      sprintf(outName,"%s/ReconstructedD0.root",path);
      TFile* outroot = new TFile(outName,"recreate");
      h1->Write();
      h3->Write();
      outroot->Close();
      delete outroot;
    }
  }
  if(PtD0){
    Char_t outName[1024];
    sprintf(outName,"%s/ReconstructedD0.root",path);
    TFile* outroot = new TFile(outName,"recreate");
    h1->Write();
    h2->Write();
    h3->Write();
    h4->Write();
    outroot->Close();
    delete outroot;
  }
  Char_t foutName[1024];
  sprintf(foutName,"%s/Cuts",path);
  ofstream fout(foutName);
  Char_t st2[1024];
  sprintf(st2,"First Pt:       %g",kPtCut);
  fout<<st2<<endl;
  sprintf(st2,"d0 low:         %g",kd0Cut);
  fout<<st2<<endl;
  sprintf(st2,"d0 high:        %g",kd0CutHigh);
  fout<<st2<<endl;
  sprintf(st2,"V0 low:         %g",cuts[0]);
  fout<<st2<<endl;
  sprintf(st2,"V0 high:        %g",cuts[1]);
  fout<<st2<<endl;
  sprintf(st2,"InvMass Diff:   %g",cuts[2]);
  fout<<st2<<endl;
  sprintf(st2,"cosPointAng:    %g",cuts[3]);
  fout<<st2<<endl;
  sprintf(st2,"d0d0:           %g",cuts[4]);
  fout<<st2<<endl;
  sprintf(st2,"cosTheta*:      %g",cuts[5]);
  fout<<st2<<endl;
  sprintf(st2,"PtChild:        %g",cuts[6]);
  fout<<st2<<endl;
  sprintf(st2,"DCA:            %g",cutDCA);
  fout<<st2<<endl;
  fout.close();

  Char_t fName[1024];
  sprintf(fName,"%s/Events",path);
  ofstream fevent(fName);
  for(int i=0;i<nD0;i++){fevent<<event[i]<<endl;}
  fevent.close();
}
//___________________________________________________________________________
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(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
  //if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;

  return kTRUE;
}
//____________________________________________________________________________
void GetPrimaryVertex(int i,Char_t* path="./") {

  int event=i;

  Char_t falice[1024];
  sprintf(falice,"%s/galice.root",path);
  TFile * galice = new TFile(falice);
  
  TDirectory * curdir;  

  Char_t vname[20];
  galice->cd();
  
  sprintf(vname,"Vertex_%d",event);
  TArrayF o = 0;
  o.Set(3);
  AliHeader * header = 0;
  
  TTree * treeE = (TTree*)gDirectory->Get("TE");
  treeE->SetBranchAddress("Header",&header);
  treeE->GetEntry(event);
  AliGenEventHeader* genHeader = header->GenEventHeader();
  if(genHeader){
    genHeader->PrimaryVertex(o);
    primaryvertex[0] = (Double_t)o[0];
    primaryvertex[1] = (Double_t)o[1];
    primaryvertex[2] = (Double_t)o[2];
  }
  else{
    printf("Can't find Header");
  }
  delete header;
  delete galice;
}
//____________________________________________________________________________
PtD0(Char_t* path="./",TH1F * h1,TH1F * h4){
  
  Char_t falice[1024];
  sprintf(falice,"%s/galice.root",path);  
  TFile *f = new TFile(falice);
  gAlice=(AliRun*)f->Get("gAlice");
  
  TParticle *p;
  int nd0=0;
  int nkminus =0;
  int npipluss = 0;
  int nEvent=gAlice->GetEventsPerRun();
  
  for (Int_t i = 0; i <nEvent; i++) {
    cout<<"Event: "<<i<<endl;
    gAlice->GetEvent(i);
    Int_t nPart = gAlice->GetNtrack();
    for (Int_t iPart = 0; iPart < nPart; iPart++) {
      //cout<<"Particlenr.: "<<iPart<<endl;
      p = (TParticle*)gAlice->Particle(iPart);
      if (p->GetPdgCode()==421){
	if(fabs(p->Eta())<0.9) h1->Fill(p.Pt());
	h4->Fill(p.Eta());
	nd0++;
      }
      if (p->GetPdgCode()==-321){
	nkminus++;
      }
      if (p->GetPdgCode()==211){
	npipluss++;
      }
    }
  }
  delete gAlice;
}
Commit	Line	Data
0bd0c1ef	1	//#define __COMPILE__
	2	//#ifdef __COMPILE__
	3	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	4	//-- --- standard headers-------------
	5	#include <iostream.h>
68772666	6	#include <fstream.h>
0bd0c1ef	7	//--------Root headers ---------------
	8	#include <TSystem.h>
	9	#include <TFile.h>
	10	#include <TString.h>
	11	#include <TStopwatch.h>
	12	#include <TObject.h>
	13	#include <TVector3.h>
	14	#include <TTree.h>
68772666	15	#include <TObjArray>
0bd0c1ef	16	#include <TParticle.h>
	17	#include <TArray.h>
	18	//----- AliRoot headers ---------------
	19	#include "alles.h"
	20	#include "AliRun.h"
	21	#include "AliKalmanTrack.h"
	22	#include "AliITStrackV2.h"
	23	#include "AliHeader.h"
	24	#include "AliGenEventHeader.h"
	25	#include "AliV0vertex.h"
	26	#include "AliV0vertexer.h"
	27	#include "AliITSVertex.h"
	28	#include "AliITSVertexer.h"
	29	#include "AliITSVertexerTracks.h"
	30	#include "AliD0Trigger.h"
	31	#endif
	32	//-------------------------------------
68772666	33	typedef struct {
	34	Int_t lab;
	35	Int_t pdg;
	36	Int_t mumlab;
	37	Int_t mumpdg;
	38	Float_t Vx,Vy,Vz;
	39	Float_t Px,Py,Pz;
	40	} RECTRACK;
	41
0bd0c1ef	42	// field (T)
	43	const Double_t kBz = 0.4;
	44
	45	// primary vertex
de3c3890	46	double primaryvertex[3]={0.,0.,0.};
0bd0c1ef	47
	48	//sec. vertex
	49	double v2[3]={0,0,0};
	50
	51	// sigle track cuts
de3c3890	52	const Double_t kPtCut = 0.5; // 0.5 GeV/c
	53	const Double_t kd0Cut = 50.; // 50 micron
	54	const Double_t kd0CutHigh = 400.; // 200 micron
0bd0c1ef	55
	56
	57	//cuts for combined tracks
de3c3890	58	const Double_t cuts[7] = {0.005, // 0.005 cuts[0] = lowest V0 cut (cm)
	59	0.800, // 0.015 cuts[1] = highest V0 cut (cm)
	60	0.012, // 0.012 cuts[2] = inv. mass cut (diferense) (Gev/c)
	61	0.8, // 0.8 cuts[3] = min. cosine value for pointing angle
	62	-5000, // -5000 cuts[4] = d0d0
68772666	63	0, // 0.8 cuts[5] = costhetastar
de3c3890	64	0.5}; // 0.5 cuts[6] = ptchild
0bd0c1ef	65	//cut for distance of closest aprach
de3c3890	66	double cutDCA=0.05; //0.05
0bd0c1ef	67
	68	// this function applies single track cuts
	69	Bool_t TrkCuts(const AliITStrackV2& trk);
	70
	71	// this function creates TObjArrays with positive and negative tracks
	72	void SelectTracks(TTree& itsTree,
	73	TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
	74	TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
	75
	76	//void GetPrimaryVertex(int i,Char_t* path="./");
	77
de3c3890	78	//void PtD0(Char_t* path="./");
de3c3890	79
0bd0c1ef	80	Int_t iTrkP,iTrkN,itsEntries;
68772666	81	Char_t trksName[100],refsName[100];
0bd0c1ef	82	Int_t nTrksP=0,nTrksN=0;
	83	Int_t nD0=0;
	84	int ev=0;
	85	double mom[6];
de3c3890	86	int event[10000];
de3c3890	87	int index=0;
68772666	88	Bool_t isSignal;
	89	Int_t nTotEv=0,nD0recSgn=0,nD0recBkgS=0,nD0recBkg=0,nD0rec1ev=0;
	90	Int_t pdg[2],mum[2],mumlab[2];
	91	RECTRACK rectrk;
0bd0c1ef	92
de3c3890	93	void RunD0offline(Char_t* path="./",bool h=false,bool PtD0=false) {
68772666	94
0bd0c1ef	95	AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
68772666	96
0bd0c1ef	97	// Open file with ITS tracks
0bd0c1ef	98	Char_t fnameTrack[1024];
de3c3890	99	//sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
0bd0c1ef	100	sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
0bd0c1ef	101	TFile* itstrks = TFile::Open(fnameTrack);
68772666	102
	103	Char_t refFile[1024];
	104	//sprintf(fnameTrack,"%s/ITStracksRefFile.root",path);
	105	sprintf(refFile,"%s/ITStracksRefFile.root",path);
	106	TFile* itsrefs = TFile::Open(refFile);
	107
0bd0c1ef	108
0bd0c1ef	109	//the offline stuff
	110	// define the cuts for vertexing
	111	Double_t vtxcuts[]={33., // max. allowed chi2
	112	0.0, // min. allowed negative daughter's impact param
	113	0.0, // min. allowed positive daughter's impact param
	114	1.0, // max. allowed DCA between the daughter tracks
de3c3890	115	-1.0, // min. allowed cosine of V0's pointing angle
0bd0c1ef	116	0.0, // min. radius of the fiducial volume
0bd0c1ef	117	2.9};// max. radius of the fiducial volume
68772666	118
de3c3890	119	TH1F *h1 = new TH1F("h1","Transvers momentun of reconstructed D0",100,0,10);
	120	TH1F *h2 = new TH1F("h2","Transvers momentun of D0 with \|Eta\|<0.9",100,0,10);
	121	TH1F *h3 = new TH1F("h3","Eta reconstructed of D0",100,-5,5);
	122	TH1F *h4 = new TH1F("h4","Eta of D0",100,-5,5);
68772666	123
	124	Char_t falice[1024];
	125	sprintf(falice,"%s/galice.root",path);
	126	TFile *f = new TFile(falice);
	127	gAlice=(AliRun*)f->Get("gAlice");
	128	int nEvent=gAlice->GetEventsPerRun();
	129	//int nEvent=20;
	130	cout<<"#Events: "<<nEvent<<endl;
	131	delete gAlice;
de3c3890	132
68772666	133	const Char_t *name="AliD0offline";
	134	cerr<<'\n'<<name<<"...\n";
	135	gBenchmark->Start(name);
	136
de3c3890	137	for(ev=0;ev<nEvent;ev++) {
	138
	139	cout<<"\n Event: "<<ev<<endl;
68772666	140
de3c3890	141	// tracks from ITS
	142	sprintf(trksName,"TreeT_ITS_%d",ev);
	143	TTree itsTree=(TTree)itstrks->Get(trksName);
	144	if(!itsTree) continue;
	145	itsEntries = (Int_t)itsTree->GetEntries();
	146	printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
	147
68772666	148	// tree from reference file
	149
	150	sprintf(refsName,"Tree_Ref_%d",ev);
	151	TTree refTree=(TTree)itsrefs->Get(refsName);
	152	refTree->SetBranchAddress("rectracks",&rectrk);
	153
de3c3890	154	//Getting primary Vertex
	155	GetPrimaryVertex(ev,path);
	156
68772666	157	// count the total number of events
	158	nTotEv++;
	159
de3c3890	160	// call function which applies sigle track selection and
	161	// separetes positives and negatives
	162	TObjArray trksP(itsEntries/2);
	163	Int_t *itsEntryP = new Int_t[itsEntries];
	164	TObjArray trksN(itsEntries/2);
	165	Int_t *itsEntryN = new Int_t[itsEntries];
	166	SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
	167
	168	cout<<"#pos: "<<nTrksP<<endl;
	169	cout<<"#neg: "<<nTrksN<<endl;
	170
	171	// create the AliV0vertexer object
	172	AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
	173
	174	AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
	175
	176	double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
	177
	178	for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
	179	postrack = (AliITStrackV2*)trksP.At(iTrkP);
68772666	180
	181	// get info on tracks PDG and mothers PDG from reference file
	182	refTree->GetEvent(itsEntryP[iTrkP]);
	183	pdg[0] = rectrk.pdg;
	184	mum[0] = rectrk.mumpdg;
	185	mumlab[0] = rectrk.mumlab;
	186
de3c3890	187	for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
de3c3890	188	negtrack = (AliITStrackV2*)trksN.At(iTrkN);
68772666	189
	190	// get info on tracks PDG and mothers PDG from reference file
	191	refTree->GetEvent(itsEntryN[iTrkN]);
	192	pdg[1] = rectrk.pdg;
	193	mum[1] = rectrk.mumpdg;
	194	mumlab[1] = rectrk.mumlab;
	195
de3c3890	196	D0->SetTracks(postrack,negtrack);
68772666	197
de3c3890	198	// ----------- DCA MINIMIZATION ------------------
	199	//
	200	// find the DCA and propagate the tracks to the DCA
	201	double dca = vertexer2->PropagateToDCA(negtrack,postrack);
	202
	203	if(dca<cutDCA){
	204	// define the AliV0vertex object
	205	AliV0vertex vertex2 = new AliV0vertex(negtrack,*postrack);
	206	// get position of the vertex
	207	vertex2->GetXYZ(v2[0],v2[1],v2[2]);
	208	delete vertex2;
	209	//if(D0->FindV0offline(v2) && D0->d0d0()){
	210	if(D0->d0d0()){
	211	D0->SetV0(v2);
	212	D0->FindMomentaOffline();
	213	if(D0->FindInvMass() && D0->CosThetaStar() && D0->PointingAngle() && D0->pTchild()){
	214	nD0++;
	215	event[index]=ev; index++;
	216	if(h){
	217	h1->Fill(D0->Pt());
	218	h3->Fill(D0->Eta());
0bd0c1ef	219	}
68772666	220
	221	if(mumlab[0]==mumlab[1] && TMath::Abs(mum[0])==421 && TMath::Abs(mum[1])==421) {
	222	nD0recSgn++;
	223	}
	224	else {
	225	if(TMath::Abs(mum[0])==421 \|\| TMath::Abs(mum[0])==411 \|\|
	226	TMath::Abs(mum[1])==421 \|\| TMath::Abs(mum[1])==411) {
	227	nD0recBkgS++;
	228	} else {
	229	nD0recBkg++;
	230	}
	231	}
de3c3890	232	}
de3c3890	233	}
0bd0c1ef	234	}
0bd0c1ef	235	}
de3c3890	236	}
68772666	237
	238	//delete D0;
	239	//delete itstrks;
	240	//delete itsTree;
	241	//delete trksP;
	242	//delete itsEntryP;
	243	//delete trksN;
	244	//delete itsEntryN;
de3c3890	245	}
68772666	246
	247	cout<<"\nMy #D0: "<<nD0<<"\n"<<endl;
	248	printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
	249	printf("\n+++\n+++ Total number of D0 candidates:\n
	250	+++ Sgn: %d\n
	251	+++ BkgS: %d\n
	252	+++ Bkg: %d\n+++\n",nD0recSgn,nD0recBkgS,nD0recBkg);
	253
0bd0c1ef	254	gBenchmark->Stop(name);
0bd0c1ef	255	gBenchmark->Show(name);
68772666	256
de3c3890	257	if(h){
	258	if(!PtD0){
	259	TCanvas *c = new TCanvas("c","c",700,1000);
	260	c->Divide(1,2);
	261	c->cd(1);
	262	h1->Draw();
	263	pt = new TPaveText(7.3,0.4,11,1.8,"br");
	264	pt->SetFillColor(0);
	265	pt->SetBorderSize(1);
	266	pt->AddText("Cuts:");
	267	Char_t st[1024];
	268	sprintf(st,"First Pt: %g",kPtCut);
	269	pt->AddText(st);
	270	sprintf(st,"d0 low: %g",kd0Cut);
	271	pt->AddText(st);
	272	sprintf(st,"d0 high: %g",kd0CutHigh);
	273	pt->AddText(st);
	274	sprintf(st,"V0 low: %g",cuts[0]);
	275	pt->AddText(st);
	276	sprintf(st,"V0 high: %g",cuts[1]);
	277	pt->AddText(st);
	278	sprintf(st,"InvMass Diff: %g",cuts[2]);
	279	pt->AddText(st);
	280	sprintf(st,"cosPointingAngle: %g",cuts[3]);
	281	pt->AddText(st);
	282	sprintf(st,"d0d0: %g",cuts[4]);
	283	pt->AddText(st);
	284	sprintf(st,"cosThetaStar: %g",cuts[5]);
	285	pt->AddText(st);
	286	sprintf(st,"PtChild: %g",cuts[6]);
	287	pt->AddText(st);
	288	sprintf(st,"DCA: %g",cutDCA);
	289	pt->AddText(st);
	290	pt->Draw();
	291	c_1->Modified();
	292	c->cd();
	293	c->cd(2);
	294	h3->Draw();
	295	}
	296	if(PtD0){
	297	PtD0(path,h2,h4);
	298	TCanvas *c = new TCanvas("c","c",1000,700);
	299	c->Divide(2,2);
	300	c->cd(1);
	301	h1->Draw();
	302	pt = new TPaveText(7.3,0.4,11,1.8,"br");
	303	pt->SetFillColor(0);
	304	pt->SetBorderSize(1);
	305	pt->AddText("Cuts:");
	306	Char_t st[1024];
	307	sprintf(st,"First Pt: %g",kPtCut);
	308	pt->AddText(st);
	309	sprintf(st,"d0 low: %g",kd0Cut);
	310	pt->AddText(st);
	311	sprintf(st,"d0 high: %g",kd0CutHigh);
	312	pt->AddText(st);
	313	sprintf(st,"V0 low: %g",cuts[0]);
	314	pt->AddText(st);
	315	sprintf(st,"V0 high: %g",cuts[1]);
	316	pt->AddText(st);
	317	sprintf(st,"InvMass Diff: %g",cuts[2]);
	318	pt->AddText(st);
68772666	319	sprintf(st,"cosPointAng: %g",cuts[3]);
de3c3890	320	pt->AddText(st);
	321	sprintf(st,"d0d0: %g",cuts[4]);
	322	pt->AddText(st);
68772666	323	sprintf(st,"cosTheta*: %g",cuts[5]);
de3c3890	324	pt->AddText(st);
	325	sprintf(st,"PtChild: %g",cuts[6]);
	326	pt->AddText(st);
	327	sprintf(st,"DCA: %g",cutDCA);
	328	pt->AddText(st);
	329	pt->Draw();
	330	c_1->Modified();
	331	c->cd();
	332	c->cd(2);
	333	h3->Draw();
	334	c->cd(3);
	335	h2->Draw();
	336	c->cd(4);
	337	h4->Draw();
68772666	338	}
	339	}
	340	if(h){
	341	if(!PtD0){
	342	Char_t outName[1024];
	343	sprintf(outName,"%s/ReconstructedD0.root",path);
	344	TFile* outroot = new TFile(outName,"recreate");
	345	h1->Write();
	346	h3->Write();
	347	outroot->Close();
	348	delete outroot;
de3c3890	349	}
68772666	350	}
	351	if(PtD0){
	352	Char_t outName[1024];
	353	sprintf(outName,"%s/ReconstructedD0.root",path);
	354	TFile* outroot = new TFile(outName,"recreate");
	355	h1->Write();
	356	h2->Write();
	357	h3->Write();
	358	h4->Write();
	359	outroot->Close();
	360	delete outroot;
	361	}
	362	Char_t foutName[1024];
	363	sprintf(foutName,"%s/Cuts",path);
	364	ofstream fout(foutName);
	365	Char_t st2[1024];
	366	sprintf(st2,"First Pt: %g",kPtCut);
	367	fout<<st2<<endl;
	368	sprintf(st2,"d0 low: %g",kd0Cut);
	369	fout<<st2<<endl;
	370	sprintf(st2,"d0 high: %g",kd0CutHigh);
	371	fout<<st2<<endl;
	372	sprintf(st2,"V0 low: %g",cuts[0]);
	373	fout<<st2<<endl;
	374	sprintf(st2,"V0 high: %g",cuts[1]);
	375	fout<<st2<<endl;
	376	sprintf(st2,"InvMass Diff: %g",cuts[2]);
	377	fout<<st2<<endl;
	378	sprintf(st2,"cosPointAng: %g",cuts[3]);
	379	fout<<st2<<endl;
	380	sprintf(st2,"d0d0: %g",cuts[4]);
	381	fout<<st2<<endl;
	382	sprintf(st2,"cosTheta*: %g",cuts[5]);
	383	fout<<st2<<endl;
	384	sprintf(st2,"PtChild: %g",cuts[6]);
	385	fout<<st2<<endl;
	386	sprintf(st2,"DCA: %g",cutDCA);
	387	fout<<st2<<endl;
	388	fout.close();
	389
	390	Char_t fName[1024];
	391	sprintf(fName,"%s/Events",path);
	392	ofstream fevent(fName);
	393	for(int i=0;i<nD0;i++){fevent<<event[i]<<endl;}
	394	fevent.close();
0bd0c1ef	395	}
	396	//___________________________________________________________________________
	397	void SelectTracks(TTree& itsTree,
	398	TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
	399	TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
	400	//
	401	// this function creates two TObjArrays with positive and negative tracks
	402	//
	403	nTrksP=0,nTrksN=0;
	404
	405
	406	Int_t entr = (Int_t)itsTree.GetEntries();
	407
	408	// trasfer tracks from tree to arrays
	409	for(Int_t i=0; i<entr; i++) {
	410
	411	AliITStrackV2 *itstrack = new AliITStrackV2;
	412	itsTree.SetBranchAddress("tracks",&itstrack);
	413
	414	itsTree.GetEvent(i);
	415
	416	// single track selection
	417	if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
	418
	419	if(itstrack->Get1Pt()>0.) { // negative track
	420	trksN.AddLast(itstrack);
	421	itsEntryN[nTrksN] = i;
	422	nTrksN++;
	423	} else { // positive track
	424	trksP.AddLast(itstrack);
	425	itsEntryP[nTrksP] = i;
	426	nTrksP++;
	427	}
	428
	429	}
	430
	431	return;
	432	}
	433	//____________________________________________________________________________
	434	Bool_t TrkCuts(const AliITStrackV2& trk) {
	435	//
	436	// this function tells if track passes some kinematical cuts
	437	//
	438	if(TMath::Abs(1./trk.Get1Pt()) < kPtCut) return kFALSE;
	439	if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
de3c3890	440	//if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
0bd0c1ef	441
	442	return kTRUE;
	443	}
	444	//____________________________________________________________________________
	445	void GetPrimaryVertex(int i,Char_t* path="./") {
	446
	447	int event=i;
	448
	449	Char_t falice[1024];
	450	sprintf(falice,"%s/galice.root",path);
	451	TFile * galice = new TFile(falice);
	452
	453	TDirectory * curdir;
	454
	455	Char_t vname[20];
	456	galice->cd();
	457
	458	sprintf(vname,"Vertex_%d",event);
	459	TArrayF o = 0;
	460	o.Set(3);
	461	AliHeader * header = 0;
	462
	463	TTree * treeE = (TTree*)gDirectory->Get("TE");
	464	treeE->SetBranchAddress("Header",&header);
	465	treeE->GetEntry(event);
	466	AliGenEventHeader* genHeader = header->GenEventHeader();
	467	if(genHeader){
	468	genHeader->PrimaryVertex(o);
	469	primaryvertex[0] = (Double_t)o[0];
	470	primaryvertex[1] = (Double_t)o[1];
	471	primaryvertex[2] = (Double_t)o[2];
	472	}
	473	else{
	474	printf("Can't find Header");
	475	}
	476	delete header;
	477	delete galice;
	478	}
de3c3890	479	//____________________________________________________________________________
	480	PtD0(Char_t* path="./",TH1F * h1,TH1F * h4){
	481
	482	Char_t falice[1024];
	483	sprintf(falice,"%s/galice.root",path);
	484	TFile *f = new TFile(falice);
	485	gAlice=(AliRun*)f->Get("gAlice");
	486
	487	TParticle *p;
	488	int nd0=0;
	489	int nkminus =0;
	490	int npipluss = 0;
	491	int nEvent=gAlice->GetEventsPerRun();
68772666	492
de3c3890	493	for (Int_t i = 0; i <nEvent; i++) {
	494	cout<<"Event: "<<i<<endl;
	495	gAlice->GetEvent(i);
	496	Int_t nPart = gAlice->GetNtrack();
	497	for (Int_t iPart = 0; iPart < nPart; iPart++) {
	498	//cout<<"Particlenr.: "<<iPart<<endl;
	499	p = (TParticle*)gAlice->Particle(iPart);
	500	if (p->GetPdgCode()==421){
	501	if(fabs(p->Eta())<0.9) h1->Fill(p.Pt());
	502	h4->Fill(p.Eta());
	503	nd0++;
	504	}
	505	if (p->GetPdgCode()==-321){
	506	nkminus++;
	507	}
	508	if (p->GetPdgCode()==211){
	509	npipluss++;
	510	}
	511	}
	512	}
	513	delete gAlice;
	514	}