[u/mrichter/AliRoot.git] / ITS / DoVerticesSPD.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include "TFile.h"
#include "TGeoManager.h"
#include "TInterpreter.h"
#include "TClassTable.h"
#include "TNtuple.h"
#include "TTree.h"
#include "TArrayF.h"
#include "TParticle.h"
#include "AliESD.h"
#include "AliGenPythiaEventHeader.h"
#include "AliTracker.h"
#include "AliHeader.h"
#include "AliITSLoader.h"
#include "AliVertexerTracks.h"
#include "AliCDBManager.h"
#include "AliGeomManager.h"
#include "AliITSVertexer3D.h"
#include "AliITSVertexerZ.h"
#include "AliESDVertex.h"
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliGenEventHeader.h"
#include "AliStack.h"
#endif

/*  $Id$    */

Bool_t DoVerticesSPD(Int_t optdebug=1){
  TFile *fint = new TFile("VertexSPD.root","recreate");
  TNtuple *nt = new TNtuple("ntvert","vertices","xtrue:ytrue:ztrue:zZ:zdiffZ:zerrZ:ntrksZ:x3D:xdiff3D:xerr3D:y3D:ydiff3D:yerr3D:z3D:zdiff3D:zerr3D:ntrks3D:dndy:ntrklets:nrp1:ptyp");

  if (gClassTable->GetID("AliRun") < 0) {
    gInterpreter->ExecuteMacro("loadlibs.C");
  }
  // Set OCDB if needed
  AliCDBManager* man = AliCDBManager::Instance();
  if (!man->IsDefaultStorageSet()) {
    printf("Setting a local default storage and run number 0\n");
    man->SetDefaultStorage("local://$ALICE_ROOT");
    man->SetRun(0);
  }
  else {
    printf("Using deafult storage \n");
  }
 
  // retrives geometry 
  if(!gGeoManager){
    AliGeomManager::LoadGeometry("geometry.root");
  }

  AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
  if (!runLoader) {
    Error("DoVertices", "getting run loader from file %s failed", 
	  "galice.root");
    return kFALSE;
  }
  runLoader->LoadgAlice();
  gAlice = runLoader->GetAliRun();
  if (!gAlice) {
    Error("DoVertices", "no galice object found");
    return kFALSE;
  }
  runLoader->LoadKinematics();
  runLoader->LoadHeader();
  AliITSLoader* ITSloader =  (AliITSLoader*) runLoader->GetLoader("ITSLoader");
  ITSloader->LoadRecPoints("read");

  AliMagF * magf = gAlice->Field();
  AliTracker::SetFieldMap(magf,kTRUE);
  if(optdebug) printf("MagneticField=%f\n",AliTracker::GetBz());
  
  Int_t totev=runLoader->GetNumberOfEvents();
  if(optdebug)  printf("Number of events= %d\n",totev);

  
//   TFile* esdFile = TFile::Open("AliESDs.root");
//   if (!esdFile || !esdFile->IsOpen()) {
//     Error("DoVertices", "opening ESD file %s failed", "AliESDs.root");
//     return kFALSE;
//   }
//   AliESD* esd = new AliESD;
//   TTree* tree = (TTree*) esdFile->Get("esdTree");
//   if (!tree) {
//     Error("DoVertices", "no ESD tree found");
//     return kFALSE;
//   }
//   tree->SetBranchAddress("ESD", &esd);

  Double_t xnom=0.,ynom=0.;
  AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
  vertz->Init("default");
  AliITSVertexer3D *vert3d = new AliITSVertexer3D();
  vert3d->Init("default");
  //  vert3d->SetDebug(10);
  //  vertz->ConfigIterations(5);

  Int_t goodz=0,good3d=0;

  for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
    TArrayF mcVertex(3); 
    runLoader->GetEvent(iEvent);
    runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
    AliGenPythiaEventHeader *evh=(AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
    Int_t ptype = evh->ProcessType();
    if(optdebug){
      printf("==============================================================\n");
      printf("\nEvent: %d ---- Process Type = %d \n",iEvent,ptype);
    }

    AliStack* stack = runLoader->Stack();
    TTree *treek=(TTree*)runLoader->TreeK();
    Int_t npart = (Int_t)treek->GetEntries();
    if(optdebug) printf("particles  %d\n",npart);

    Double_t dNchdy = 0.;

   // loop on particles
    for(Int_t pa=0; pa<npart; pa++) {
      TParticle* part = (TParticle*)stack->Particle(pa);
      Int_t pdg = part->GetPdgCode();
      Int_t apdg = TMath::Abs(pdg);
      Double_t energy  = part->Energy();
      if(energy>6900.) continue; // reject incoming protons
      Double_t pz = part->Pz();
      Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));


      if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;      // reject secondaries
      if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
      if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
    }
    if(optdebug) printf(" dNch/dy = %f\n",dNchdy);
 
    TTree* cltree = ITSloader->TreeR();

    AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
    AliMultiplicity *alimult = vertz->GetMultiplicity();
    Int_t ntrklets=0,nrecp1=0;
    if(alimult) {
      nrecp1=alimult->GetNumberOfTracklets() ;
      ntrklets = alimult->GetNumberOfTracklets();
      for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
	if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
      }
    }

    AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);

    TDirectory *current = gDirectory;
    fint->cd();
    Float_t xnt[21];

    Double_t zz = 0.;
    Double_t zresz = 0.;
    Double_t zdiffz = 0.;

    Int_t ntrkz = 0;
    if(vtxz){
 
      ntrkz = vtxz->GetNContributors();
      if(ntrkz>0)goodz++;
      zz=vtxz->GetZv();
      zresz =vtxz->GetZRes(); // microns
      zdiffz = 10000.*(zz-mcVertex[2]); // microns
    }

    Double_t x3d = 0.;
    Double_t xerr3d = 0.;
    Double_t xdiff3d = 0.;

    Double_t y3d = 0.;
    Double_t yerr3d = 0.;
    Double_t ydiff3d = 0.;

    Double_t z3d = 0.;
    Double_t zerr3d = 0.;
    Double_t zdiff3d = 0.;

    Int_t ntrk3d = 0;
    if(vtx3d){

      ntrk3d = vtx3d->GetNContributors();
      if(ntrk3d>0)good3d++;
      x3d = vtx3d->GetXv();
      xerr3d=vtx3d->GetXRes();
      xdiff3d = 10000.*(x3d-mcVertex[0]);  // microns

      y3d = vtx3d->GetYv();
      yerr3d=vtx3d->GetYRes();
      ydiff3d = 10000.*(y3d-mcVertex[1]);  // microns

      z3d = vtx3d->GetZv();
      zerr3d=vtx3d->GetZRes();
      zdiff3d = 10000.*(z3d-mcVertex[2]);  // microns

    }

    xnt[0]=mcVertex[0];//x
    xnt[1]=mcVertex[1];//y
    xnt[2]=mcVertex[2];//z

    xnt[3]=zz;
    xnt[4]=zdiffz;
    xnt[5]=zresz;
    xnt[6]=ntrkz;

    xnt[7]=x3d;
    xnt[8]=xdiff3d;
    xnt[9]=xerr3d;
    xnt[11]=y3d;
    xnt[11]=ydiff3d;
    xnt[12]=yerr3d;
    xnt[13]=z3d;
    xnt[14]=zdiff3d;
    xnt[15]=zerr3d;
    xnt[16]=ntrk3d;

    xnt[17]=dNchdy;
    xnt[18]=float(ntrklets);
    xnt[19]=float(nrecp1);
    xnt[20]=float(ptype);
    nt->Fill(xnt);
    current->cd();
    
    if(optdebug){
      printf("\nVertexerZ:  \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",zz,ntrkz,mcVertex[2],zdiffz);
      printf("Vertexer3D: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",z3d,ntrk3d,mcVertex[2],zdiff3d);
    }
    
  }
  fint->cd();
  nt->Write();
  fint->Close();
  delete fint;
  if(optdebug){
   printf("***********************************************\n");
   printf("Number of Z vertices found= %d\n",goodz);
   printf("efficiency=%f\n",float(goodz)/float(totev));
   printf("Number of 3D vertices found= %d\n",good3d);
   printf("efficiency=%f\n",float(good3d)/float(totev));
  }
  return kTRUE;
}
Commit	Line	Data
f5826889	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include "TFile.h"
	3	#include "TGeoManager.h"
	4	#include "TInterpreter.h"
	5	#include "TClassTable.h"
	6	#include "TNtuple.h"
	7	#include "TTree.h"
	8	#include "TArrayF.h"
	9	#include "TParticle.h"
	10	#include "AliESD.h"
	11	#include "AliGenPythiaEventHeader.h"
	12	#include "AliTracker.h"
	13	#include "AliHeader.h"
	14	#include "AliITSLoader.h"
	15	#include "AliVertexerTracks.h"
	16	#include "AliCDBManager.h"
	17	#include "AliGeomManager.h"
	18	#include "AliITSVertexer3D.h"
	19	#include "AliITSVertexerZ.h"
	20	#include "AliESDVertex.h"
	21	#include "AliRun.h"
	22	#include "AliRunLoader.h"
	23	#include "AliGenEventHeader.h"
	24	#include "AliStack.h"
	25	#endif
	26
	27	/* $Id$ */
	28
	29	Bool_t DoVerticesSPD(Int_t optdebug=1){
	30	TFile *fint = new TFile("VertexSPD.root","recreate");
	31	TNtuple *nt = new TNtuple("ntvert","vertices","xtrue:ytrue:ztrue:zZ:zdiffZ:zerrZ:ntrksZ:x3D:xdiff3D:xerr3D:y3D:ydiff3D:yerr3D:z3D:zdiff3D:zerr3D:ntrks3D:dndy:ntrklets:nrp1:ptyp");
	32
	33	if (gClassTable->GetID("AliRun") < 0) {
	34	gInterpreter->ExecuteMacro("loadlibs.C");
	35	}
	36	// Set OCDB if needed
	37	AliCDBManager* man = AliCDBManager::Instance();
	38	if (!man->IsDefaultStorageSet()) {
	39	printf("Setting a local default storage and run number 0\n");
	40	man->SetDefaultStorage("local://$ALICE_ROOT");
	41	man->SetRun(0);
	42	}
	43	else {
	44	printf("Using deafult storage \n");
	45	}
	46
	47	// retrives geometry
	48	if(!gGeoManager){
	49	AliGeomManager::LoadGeometry("geometry.root");
	50	}
	51
	52	AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
	53	if (!runLoader) {
	54	Error("DoVertices", "getting run loader from file %s failed",
	55	"galice.root");
	56	return kFALSE;
	57	}
	58	runLoader->LoadgAlice();
	59	gAlice = runLoader->GetAliRun();
	60	if (!gAlice) {
	61	Error("DoVertices", "no galice object found");
	62	return kFALSE;
	63	}
	64	runLoader->LoadKinematics();
65	runLoader->LoadHeader();
66	AliITSLoader* ITSloader = (AliITSLoader*) runLoader->GetLoader("ITSLoader");
67	ITSloader->LoadRecPoints("read");
68
69	AliMagF * magf = gAlice->Field();
70	AliTracker::SetFieldMap(magf,kTRUE);
71	if(optdebug) printf("MagneticField=%f\n",AliTracker::GetBz());
72
73	Int_t totev=runLoader->GetNumberOfEvents();
74	if(optdebug) printf("Number of events= %d\n",totev);
75
76
77	// TFile* esdFile = TFile::Open("AliESDs.root");
78	// if (!esdFile \|\| !esdFile->IsOpen()) {
79	// Error("DoVertices", "opening ESD file %s failed", "AliESDs.root");
80	// return kFALSE;
81	// }
82	// AliESD* esd = new AliESD;
83	// TTree* tree = (TTree*) esdFile->Get("esdTree");
84	// if (!tree) {
85	// Error("DoVertices", "no ESD tree found");
86	// return kFALSE;
87	// }
88	// tree->SetBranchAddress("ESD", &esd);
89
90	Double_t xnom=0.,ynom=0.;
308c2f7c	91	AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
	92	vertz->Init("default");
	93	AliITSVertexer3D *vert3d = new AliITSVertexer3D();
	94	vert3d->Init("default");
f5826889	95	// vert3d->SetDebug(10);
	96	// vertz->ConfigIterations(5);
	97
	98	Int_t goodz=0,good3d=0;
	99
	100	for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
	101	TArrayF mcVertex(3);
	102	runLoader->GetEvent(iEvent);
	103	runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
	104	AliGenPythiaEventHeader evh=(AliGenPythiaEventHeader)runLoader->GetHeader()->GenEventHeader();
	105	Int_t ptype = evh->ProcessType();
	106	if(optdebug){
	107	printf("==============================================================\n");
	108	printf("\nEvent: %d ---- Process Type = %d \n",iEvent,ptype);
	109	}
	110
	111	AliStack* stack = runLoader->Stack();
	112	TTree treek=(TTree)runLoader->TreeK();
	113	Int_t npart = (Int_t)treek->GetEntries();
	114	if(optdebug) printf("particles %d\n",npart);
	115
	116	Double_t dNchdy = 0.;
	117
	118	// loop on particles
	119	for(Int_t pa=0; pa<npart; pa++) {
	120	TParticle* part = (TParticle*)stack->Particle(pa);
	121	Int_t pdg = part->GetPdgCode();
	122	Int_t apdg = TMath::Abs(pdg);
	123	Double_t energy = part->Energy();
	124	if(energy>6900.) continue; // reject incoming protons
	125	Double_t pz = part->Pz();
	126	Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
	127
	128
	129	if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue; // reject secondaries
	130	if(TMath::Sqrt((part->Vx()-mcVertex[0])(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])(part->Vy()-mcVertex[1]))>0.0010) continue;
	131	if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in \|y\|<1
	132	}
	133	if(optdebug) printf(" dNch/dy = %f\n",dNchdy);
	134
308c2f7c	135	TTree* cltree = ITSloader->TreeR();
	136
	137	AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
f5826889	138	AliMultiplicity *alimult = vertz->GetMultiplicity();
	139	Int_t ntrklets=0,nrecp1=0;
	140	if(alimult) {
	141	nrecp1=alimult->GetNumberOfTracklets() ;
	142	ntrklets = alimult->GetNumberOfTracklets();
	143	for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
	144	if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
	145	}
	146	}
	147
308c2f7c	148	AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
f5826889	149
	150	TDirectory *current = gDirectory;
	151	fint->cd();
	152	Float_t xnt[21];
	153
	154	Double_t zz = 0.;
	155	Double_t zresz = 0.;
	156	Double_t zdiffz = 0.;
	157
	158	Int_t ntrkz = 0;
	159	if(vtxz){
	160
	161	ntrkz = vtxz->GetNContributors();
	162	if(ntrkz>0)goodz++;
	163	zz=vtxz->GetZv();
	164	zresz =vtxz->GetZRes(); // microns
	165	zdiffz = 10000.*(zz-mcVertex[2]); // microns
	166	}
	167
	168	Double_t x3d = 0.;
	169	Double_t xerr3d = 0.;
	170	Double_t xdiff3d = 0.;
	171
	172	Double_t y3d = 0.;
	173	Double_t yerr3d = 0.;
	174	Double_t ydiff3d = 0.;
	175
	176	Double_t z3d = 0.;
	177	Double_t zerr3d = 0.;
	178	Double_t zdiff3d = 0.;
	179
	180	Int_t ntrk3d = 0;
	181	if(vtx3d){
	182
	183	ntrk3d = vtx3d->GetNContributors();
	184	if(ntrk3d>0)good3d++;
	185	x3d = vtx3d->GetXv();
	186	xerr3d=vtx3d->GetXRes();
	187	xdiff3d = 10000.*(x3d-mcVertex[0]); // microns
	188
	189	y3d = vtx3d->GetYv();
	190	yerr3d=vtx3d->GetYRes();
	191	ydiff3d = 10000.*(y3d-mcVertex[1]); // microns
	192
	193	z3d = vtx3d->GetZv();
	194	zerr3d=vtx3d->GetZRes();
	195	zdiff3d = 10000.*(z3d-mcVertex[2]); // microns
	196
	197	}
	198
	199	xnt[0]=mcVertex[0];//x
	200	xnt[1]=mcVertex[1];//y
	201	xnt[2]=mcVertex[2];//z
	202
	203	xnt[3]=zz;
	204	xnt[4]=zdiffz;
	205	xnt[5]=zresz;
	206	xnt[6]=ntrkz;
	207
	208	xnt[7]=x3d;
	209	xnt[8]=xdiff3d;
	210	xnt[9]=xerr3d;
	211	xnt[11]=y3d;
	212	xnt[11]=ydiff3d;
213	xnt[12]=yerr3d;
214	xnt[13]=z3d;
215	xnt[14]=zdiff3d;
216	xnt[15]=zerr3d;
217	xnt[16]=ntrk3d;
218
219	xnt[17]=dNchdy;
220	xnt[18]=float(ntrklets);
221	xnt[19]=float(nrecp1);
222	xnt[20]=float(ptype);
223	nt->Fill(xnt);
224	current->cd();
225
226	if(optdebug){
227	printf("\nVertexerZ: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",zz,ntrkz,mcVertex[2],zdiffz);
228	printf("Vertexer3D: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",z3d,ntrk3d,mcVertex[2],zdiff3d);
229	}
230
231	}
232	fint->cd();
233	nt->Write();
234	fint->Close();
235	delete fint;
236	if(optdebug){
237	printf("***********************************************\n");
238	printf("Number of Z vertices found= %d\n",goodz);
239	printf("efficiency=%f\n",float(goodz)/float(totev));
240	printf("Number of 3D vertices found= %d\n",good3d);
241	printf("efficiency=%f\n",float(good3d)/float(totev));
242	}
243	return kTRUE;
244	}