[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 "AliITSDetTypeRec.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:is3D");

  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/OCDB");
    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);

  AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
 
  Double_t xnom=0.,ynom=0.;
  AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
  vertz->Init("default");
  AliITSVertexer3D *vert3d = new AliITSVertexer3D();
  vert3d->Init("default");
  vert3d->SetWideFiducialRegion(40.,1.);
  vert3d->PrintStatus();
  vertz->SetDetTypeRec(detTypeRec);
  vert3d->SetDetTypeRec(detTypeRec);

  /* uncomment these lines to use diamond constrain */
//   Double_t posdiam[3]={0.03,0.1,0.};
//   Double_t sigdiam[3]={0.01,0.01,10.0};
//   AliESDVertex* diam=new AliESDVertex(posdiam,sigdiam);
//   vertz->SetVtxStart(diam);
//   vert3d->SetVtxStart(diam);
  /* end lines to be uncommented to use diamond constrain */

  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 to get generated dN/dy
    for(Int_t iPart=0; iPart<npart; iPart++) {
      if(!stack->IsPhysicalPrimary(iPart)) continue;
      TParticle* part = (TParticle*)stack->Particle(iPart);
      if(part->GetPDG()->Charge() == 0) continue;
      Double_t eta=part->Eta();

      if(TMath::Abs(eta)<1.5) dNchdy+=1.; 
    }
    if(optdebug) printf(" dNch/dy = %f\n",dNchdy);
 
    TTree* cltree = ITSloader->TreeR();

    AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
    AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
    AliMultiplicity *alimult = vert3d->GetMultiplicity();
    Int_t  ntrklets=0;
    Int_t nrecp1=0;
    if(alimult){
      ntrklets=alimult->GetNumberOfTracklets() ;
      nrecp1=ntrklets+alimult->GetNumberOfSingleClusters();
    }
    TString tit=vtx3d->GetTitle();
    Bool_t is3d=kFALSE;
    if(tit.Contains("3D")) is3d=kTRUE;
       

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

    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(is3d && 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[10]=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);
    xnt[21]=float(is3d);
    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"
fdc48e5f	18	#include "AliITSDetTypeRec.h"
f5826889	19	#include "AliITSVertexer3D.h"
	20	#include "AliITSVertexerZ.h"
	21	#include "AliESDVertex.h"
	22	#include "AliRun.h"
	23	#include "AliRunLoader.h"
	24	#include "AliGenEventHeader.h"
	25	#include "AliStack.h"
	26	#endif
	27
	28	/* $Id$ */
	29
	30	Bool_t DoVerticesSPD(Int_t optdebug=1){
792b611a	31
f5826889	32	TFile *fint = new TFile("VertexSPD.root","recreate");
792b611a	33	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:is3D");
f5826889	34
	35	if (gClassTable->GetID("AliRun") < 0) {
	36	gInterpreter->ExecuteMacro("loadlibs.C");
	37	}
	38	// Set OCDB if needed
	39	AliCDBManager* man = AliCDBManager::Instance();
	40	if (!man->IsDefaultStorageSet()) {
	41	printf("Setting a local default storage and run number 0\n");
162637e4	42	man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
f5826889	43	man->SetRun(0);
	44	}
	45	else {
	46	printf("Using deafult storage \n");
	47	}
	48
	49	// retrives geometry
	50	if(!gGeoManager){
	51	AliGeomManager::LoadGeometry("geometry.root");
	52	}
	53
	54	AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
	55	if (!runLoader) {
	56	Error("DoVertices", "getting run loader from file %s failed",
	57	"galice.root");
	58	return kFALSE;
	59	}
	60	runLoader->LoadgAlice();
	61	gAlice = runLoader->GetAliRun();
	62	if (!gAlice) {
	63	Error("DoVertices", "no galice object found");
	64	return kFALSE;
	65	}
	66	runLoader->LoadKinematics();
	67	runLoader->LoadHeader();
	68	AliITSLoader* ITSloader = (AliITSLoader*) runLoader->GetLoader("ITSLoader");
	69	ITSloader->LoadRecPoints("read");
	70
	71	AliMagF * magf = gAlice->Field();
	72	AliTracker::SetFieldMap(magf,kTRUE);
	73	if(optdebug) printf("MagneticField=%f\n",AliTracker::GetBz());
	74
	75	Int_t totev=runLoader->GetNumberOfEvents();
	76	if(optdebug) printf("Number of events= %d\n",totev);
	77
fdc48e5f	78	AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
fdc48e5f	79
f5826889	80	Double_t xnom=0.,ynom=0.;
308c2f7c	81	AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
	82	vertz->Init("default");
	83	AliITSVertexer3D *vert3d = new AliITSVertexer3D();
	84	vert3d->Init("default");
792b611a	85	vert3d->SetWideFiducialRegion(40.,1.);
792b611a	86	vert3d->PrintStatus();
fdc48e5f	87	vertz->SetDetTypeRec(detTypeRec);
fdc48e5f	88	vert3d->SetDetTypeRec(detTypeRec);
ecce5370	89
	90	/* uncomment these lines to use diamond constrain */
	91	// Double_t posdiam[3]={0.03,0.1,0.};
	92	// Double_t sigdiam[3]={0.01,0.01,10.0};
	93	// AliESDVertex* diam=new AliESDVertex(posdiam,sigdiam);
	94	// vertz->SetVtxStart(diam);
	95	// vert3d->SetVtxStart(diam);
	96	/* end lines to be uncommented to use diamond constrain */
f5826889	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
ecce5370	118	// loop on particles to get generated dN/dy
	119	for(Int_t iPart=0; iPart<npart; iPart++) {
	120	if(!stack->IsPhysicalPrimary(iPart)) continue;
	121	TParticle* part = (TParticle*)stack->Particle(iPart);
	122	if(part->GetPDG()->Charge() == 0) continue;
	123	Double_t eta=part->Eta();
f5826889	124
ecce5370	125	if(TMath::Abs(eta)<1.5) dNchdy+=1.;
f5826889	126	}
	127	if(optdebug) printf(" dNch/dy = %f\n",dNchdy);
	128
308c2f7c	129	TTree* cltree = ITSloader->TreeR();
	130
	131	AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
ecce5370	132	AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
ecce5370	133	AliMultiplicity *alimult = vert3d->GetMultiplicity();
792b611a	134	Int_t ntrklets=0;
	135	Int_t nrecp1=0;
	136	if(alimult){
	137	ntrklets=alimult->GetNumberOfTracklets() ;
	138	nrecp1=ntrklets+alimult->GetNumberOfSingleClusters();
f5826889	139	}
792b611a	140	TString tit=vtx3d->GetTitle();
	141	Bool_t is3d=kFALSE;
	142	if(tit.Contains("3D")) is3d=kTRUE;
	143
ecce5370	144
f5826889	145
	146	TDirectory *current = gDirectory;
	147	fint->cd();
792b611a	148	Float_t xnt[22];
f5826889	149
	150	Double_t zz = 0.;
	151	Double_t zresz = 0.;
	152	Double_t zdiffz = 0.;
	153
	154	Int_t ntrkz = 0;
	155	if(vtxz){
	156
	157	ntrkz = vtxz->GetNContributors();
	158	if(ntrkz>0)goodz++;
	159	zz=vtxz->GetZv();
	160	zresz =vtxz->GetZRes(); // microns
	161	zdiffz = 10000.*(zz-mcVertex[2]); // microns
	162	}
	163
	164	Double_t x3d = 0.;
	165	Double_t xerr3d = 0.;
	166	Double_t xdiff3d = 0.;
	167
	168	Double_t y3d = 0.;
	169	Double_t yerr3d = 0.;
	170	Double_t ydiff3d = 0.;
	171
	172	Double_t z3d = 0.;
	173	Double_t zerr3d = 0.;
	174	Double_t zdiff3d = 0.;
	175
	176	Int_t ntrk3d = 0;
	177	if(vtx3d){
	178
	179	ntrk3d = vtx3d->GetNContributors();
792b611a	180	if(is3d && ntrk3d>0)good3d++;
f5826889	181	x3d = vtx3d->GetXv();
	182	xerr3d=vtx3d->GetXRes();
	183	xdiff3d = 10000.*(x3d-mcVertex[0]); // microns
	184
	185	y3d = vtx3d->GetYv();
	186	yerr3d=vtx3d->GetYRes();
	187	ydiff3d = 10000.*(y3d-mcVertex[1]); // microns
	188
	189	z3d = vtx3d->GetZv();
	190	zerr3d=vtx3d->GetZRes();
	191	zdiff3d = 10000.*(z3d-mcVertex[2]); // microns
	192
	193	}
	194
	195	xnt[0]=mcVertex[0];//x
	196	xnt[1]=mcVertex[1];//y
	197	xnt[2]=mcVertex[2];//z
	198
	199	xnt[3]=zz;
	200	xnt[4]=zdiffz;
	201	xnt[5]=zresz;
	202	xnt[6]=ntrkz;
	203
	204	xnt[7]=x3d;
	205	xnt[8]=xdiff3d;
	206	xnt[9]=xerr3d;
8686a665	207	xnt[10]=y3d;
f5826889	208	xnt[11]=ydiff3d;
	209	xnt[12]=yerr3d;
	210	xnt[13]=z3d;
	211	xnt[14]=zdiff3d;
	212	xnt[15]=zerr3d;
	213	xnt[16]=ntrk3d;
	214
	215	xnt[17]=dNchdy;
	216	xnt[18]=float(ntrklets);
	217	xnt[19]=float(nrecp1);
	218	xnt[20]=float(ptype);
792b611a	219	xnt[21]=float(is3d);
f5826889	220	nt->Fill(xnt);
	221	current->cd();
	222
	223	if(optdebug){
	224	printf("\nVertexerZ: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",zz,ntrkz,mcVertex[2],zdiffz);
	225	printf("Vertexer3D: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",z3d,ntrk3d,mcVertex[2],zdiff3d);
	226	}
	227
	228	}
	229	fint->cd();
	230	nt->Write();
	231	fint->Close();
	232	delete fint;
	233	if(optdebug){
	234	printf("***********************************************\n");
	235	printf("Number of Z vertices found= %d\n",goodz);
	236	printf("efficiency=%f\n",float(goodz)/float(totev));
	237	printf("Number of 3D vertices found= %d\n",good3d);
	238	printf("efficiency=%f\n",float(good3d)/float(totev));
	239	}
	240	return kTRUE;
	241	}