[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 "AliITSsegmentationSPD.h"
#include "AliVertexerTracks.h"
#include "AliCDBManager.h"
#include "AliGeomManager.h"
#include "AliGRPManager.h"
#include "AliITSDetTypeRec.h"
#include "AliITSVertexer3D.h"
#include "AliITSVertexerZ.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliGenEventHeader.h"
#include "AliStack.h"
#endif

/*  $Id$    */

Bool_t DoVerticesSPD(Bool_t isMC=kFALSE, Int_t pileupalgo=1, Int_t optdebug=0){

  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:nchgen:ntrklets:nrp1:ptyp:is3D:isTriggered");

  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->SetSpecificStorage("GRP/GRP/Data",Form("local://%s",gSystem->pwd()));
    man->SetRun(0);
  }
  else {
    printf("Using deafult storage \n");
  }
 
  // retrives geometry 
  if(!gGeoManager){
    AliGeomManager::LoadGeometry("geometry.root");
  }
  AliGeomManager::ApplyAlignObjsFromCDB("ITS");

  AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
  if (!runLoader) {
    Error("DoVertices", "getting run loader from file %s failed", 
	  "galice.root");
    return kFALSE;
  }
  
  if(isMC){
    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");

  Int_t totev=runLoader->GetNumberOfEvents();
  if(optdebug)  printf("Number of events= %d\n",totev);

  TFile* esdFile = TFile::Open("AliESDs.root");
  if (!esdFile || !esdFile->IsOpen()) {
    printf("Error in opening ESD file");
    return kFALSE;
  }

  AliESDEvent * esd = new AliESDEvent;
  TTree* tree = (TTree*) esdFile->Get("esdTree");
  if (!tree) {
    printf("Error: no ESD tree found");
    return kFALSE;
  }
  esd->ReadFromTree(tree);
  AliGRPManager * grpMan=new AliGRPManager();
  grpMan->ReadGRPEntry();
  grpMan->SetMagField();
  printf("Magnetic field set to %f T\n",AliTracker::GetBz()/10.);

  AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
  AliITSsegmentation* seg = new AliITSsegmentationSPD();  
  detTypeRec->SetSegmentationModel(0,seg);

  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->SetPileupAlgo(pileupalgo);
  vert3d->PrintStatus();
  vertz->SetDetTypeRec(detTypeRec);
  vert3d->SetDetTypeRec(detTypeRec);
  vert3d->SetComputeMultiplicity(kTRUE);
  /* 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;
  // Trigger mask
  ULong64_t spdFO = (1 << 14);
  ULong64_t v0left = (1 << 11);
  ULong64_t v0right = (1 << 12);

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

      AliStack* stack = runLoader->Stack();
      TTree *treek=(TTree*)runLoader->TreeK();
      Int_t npart = (Int_t)treek->GetEntries();
      if(optdebug) printf("Process Type = %d --- Particles  %d\n",ptype,npart);
      
      // 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) nChGen+=1.; 
      }
      if(optdebug) printf(" N. of generated charged primaries in |eta|<1.5 = %f\n",nChGen);
    }
    tree->GetEvent(iEvent);

    TTree* cltree = ITSloader->TreeR();
    ULong64_t triggerMask=esd->GetTriggerMask();
    // MB1: SPDFO || V0L || V0R
    Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)));
    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();
    }


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

    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;
    Bool_t is3d=kFALSE;
    if(vtx3d){

      if(vtx3d->IsFromVertexer3D()) is3d=kTRUE;
      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]=nChGen;
    xnt[18]=float(ntrklets);
    xnt[19]=float(nrecp1);
    xnt[20]=float(ptype);
    xnt[21]=float(is3d);
    xnt[22]=float(eventTriggered);
    nt->Fill(xnt);
    current->cd();
    
    if(optdebug){
      printf("Vertexer3D: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tContributors=%d \n",x3d,y3d,z3d,ntrk3d);
      printf("VertexerZ:  \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tContributors=%d \n",0.,0.,zz,ntrkz);
      if(isMC) printf("True Pos.: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tdN/dy=%.1f\n",mcVertex[0],mcVertex[1],mcVertex[2],nChGen);
      printf("Multiplicity: Tracklets=%d  ClustersLay1=%d\n",ntrklets,nrecp1);
    }
    
  }
  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
	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 "AliITSsegmentationSPD.h"
	16	#include "AliVertexerTracks.h"
	17	#include "AliCDBManager.h"
	18	#include "AliGeomManager.h"
	19	#include "AliGRPManager.h"
	20	#include "AliITSDetTypeRec.h"
	21	#include "AliITSVertexer3D.h"
	22	#include "AliITSVertexerZ.h"
	23	#include "AliESDVertex.h"
	24	#include "AliESDEvent.h"
	25	#include "AliRun.h"
	26	#include "AliRunLoader.h"
	27	#include "AliGenEventHeader.h"
	28	#include "AliStack.h"
	29	#endif
	30
	31	/* $Id$ */
	32
	33	Bool_t DoVerticesSPD(Bool_t isMC=kFALSE, Int_t pileupalgo=1, Int_t optdebug=0){
	34
	35	TFile *fint = new TFile("VertexSPD.root","recreate");
	36	TNtuple *nt = new TNtuple("ntvert","vertices","xtrue:ytrue:ztrue:zZ:zdiffZ:zerrZ:ntrksZ:x3D:xdiff3D:xerr3D:y3D:ydiff3D:yerr3D:z3D:zdiff3D:zerr3D:ntrks3D:nchgen:ntrklets:nrp1:ptyp:is3D:isTriggered");
	37
	38	if (gClassTable->GetID("AliRun") < 0) {
	39	gInterpreter->ExecuteMacro("loadlibs.C");
	40	}
	41	// Set OCDB if needed
	42	AliCDBManager* man = AliCDBManager::Instance();
	43	if (!man->IsDefaultStorageSet()) {
	44	printf("Setting a local default storage and run number 0\n");
	45	man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
	46	// man->SetSpecificStorage("GRP/GRP/Data",Form("local://%s",gSystem->pwd()));
	47	man->SetRun(0);
	48	}
	49	else {
	50	printf("Using deafult storage \n");
	51	}
	52
	53	// retrives geometry
	54	if(!gGeoManager){
	55	AliGeomManager::LoadGeometry("geometry.root");
	56	}
	57	AliGeomManager::ApplyAlignObjsFromCDB("ITS");
	58
	59	AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
	60	if (!runLoader) {
	61	Error("DoVertices", "getting run loader from file %s failed",
	62	"galice.root");
	63	return kFALSE;
	64	}
	65
	66	if(isMC){
	67	runLoader->LoadgAlice();
	68	gAlice = runLoader->GetAliRun();
	69	if (!gAlice) {
	70	Error("DoVertices", "no galice object found");
	71	return kFALSE;
	72	}
	73	runLoader->LoadKinematics();
	74	runLoader->LoadHeader();
	75	}
	76	AliITSLoader* ITSloader = (AliITSLoader*) runLoader->GetLoader("ITSLoader");
	77	ITSloader->LoadRecPoints("read");
	78
	79	Int_t totev=runLoader->GetNumberOfEvents();
	80	if(optdebug) printf("Number of events= %d\n",totev);
	81
	82	TFile* esdFile = TFile::Open("AliESDs.root");
	83	if (!esdFile \|\| !esdFile->IsOpen()) {
	84	printf("Error in opening ESD file");
	85	return kFALSE;
	86	}
	87
	88	AliESDEvent * esd = new AliESDEvent;
	89	TTree* tree = (TTree*) esdFile->Get("esdTree");
	90	if (!tree) {
	91	printf("Error: no ESD tree found");
	92	return kFALSE;
	93	}
	94	esd->ReadFromTree(tree);
	95	AliGRPManager * grpMan=new AliGRPManager();
	96	grpMan->ReadGRPEntry();
	97	grpMan->SetMagField();
	98	printf("Magnetic field set to %f T\n",AliTracker::GetBz()/10.);
	99
	100	AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
	101	AliITSsegmentation* seg = new AliITSsegmentationSPD();
	102	detTypeRec->SetSegmentationModel(0,seg);
	103
	104	Double_t xnom=0.,ynom=0.;
	105	AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
	106	vertz->Init("default");
	107	AliITSVertexer3D *vert3d = new AliITSVertexer3D();
	108	vert3d->Init("default");
	109	vert3d->SetWideFiducialRegion(40.,1.);
	110	vert3d->SetPileupAlgo(pileupalgo);
	111	vert3d->PrintStatus();
	112	vertz->SetDetTypeRec(detTypeRec);
	113	vert3d->SetDetTypeRec(detTypeRec);
	114	vert3d->SetComputeMultiplicity(kTRUE);
	115	/* uncomment these lines to use diamond constrain */
	116	// Double_t posdiam[3]={0.03,0.1,0.};
	117	// Double_t sigdiam[3]={0.01,0.01,10.0};
	118	// AliESDVertex* diam=new AliESDVertex(posdiam,sigdiam);
	119	// vertz->SetVtxStart(diam);
	120	// vert3d->SetVtxStart(diam);
	121	/* end lines to be uncommented to use diamond constrain */
	122
	123	Int_t goodz=0,good3d=0;
	124	// Trigger mask
	125	ULong64_t spdFO = (1 << 14);
	126	ULong64_t v0left = (1 << 11);
	127	ULong64_t v0right = (1 << 12);
	128
	129	for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
	130	TArrayF mcVertex(3);
	131	runLoader->GetEvent(iEvent);
	132	Double_t nChGen = 0.;
	133	Int_t ptype = 0;
	134	if(optdebug){
	135	printf("==============================================================\n");
	136	printf("Event: %d \n",iEvent);
	137	}
	138	if(isMC){
	139	runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
	140	AliGenPythiaEventHeader evh=(AliGenPythiaEventHeader)runLoader->GetHeader()->GenEventHeader();
	141	ptype = evh->ProcessType();
	142
	143	AliStack* stack = runLoader->Stack();
	144	TTree treek=(TTree)runLoader->TreeK();
	145	Int_t npart = (Int_t)treek->GetEntries();
	146	if(optdebug) printf("Process Type = %d --- Particles %d\n",ptype,npart);
	147
	148	// loop on particles to get generated dN/dy
	149	for(Int_t iPart=0; iPart<npart; iPart++) {
	150	if(!stack->IsPhysicalPrimary(iPart)) continue;
	151	TParticle* part = (TParticle*)stack->Particle(iPart);
	152	if(part->GetPDG()->Charge() == 0) continue;
	153	Double_t eta=part->Eta();
	154
	155	if(TMath::Abs(eta)<1.5) nChGen+=1.;
	156	}
	157	if(optdebug) printf(" N. of generated charged primaries in \|eta\|<1.5 = %f\n",nChGen);
	158	}
	159	tree->GetEvent(iEvent);
	160
	161	TTree* cltree = ITSloader->TreeR();
	162	ULong64_t triggerMask=esd->GetTriggerMask();
	163	// MB1: SPDFO \|\| V0L \|\| V0R
	164	Bool_t eventTriggered = (triggerMask & spdFO \|\| ((triggerMask & v0left) \|\| (triggerMask & v0right)));
	165	AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
	166	AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
	167	AliMultiplicity *alimult = vert3d->GetMultiplicity();
	168	Int_t ntrklets=0;
	169	Int_t nrecp1=0;
	170	if(alimult){
	171	ntrklets=alimult->GetNumberOfTracklets() ;
	172	nrecp1=ntrklets+alimult->GetNumberOfSingleClusters();
	173	}
	174
	175
	176	TDirectory *current = gDirectory;
	177	fint->cd();
	178	Float_t xnt[23];
	179
	180	Double_t zz = 0.;
	181	Double_t zresz = 0.;
	182	Double_t zdiffz = 0.;
	183
	184	Int_t ntrkz = 0;
	185	if(vtxz){
	186
	187	ntrkz = vtxz->GetNContributors();
	188	if(ntrkz>0)goodz++;
	189	zz=vtxz->GetZv();
	190	zresz =vtxz->GetZRes(); // microns
	191	zdiffz = 10000.*(zz-mcVertex[2]); // microns
	192	}
	193
	194	Double_t x3d = 0.;
	195	Double_t xerr3d = 0.;
	196	Double_t xdiff3d = 0.;
	197
	198	Double_t y3d = 0.;
	199	Double_t yerr3d = 0.;
	200	Double_t ydiff3d = 0.;
	201
	202	Double_t z3d = 0.;
	203	Double_t zerr3d = 0.;
	204	Double_t zdiff3d = 0.;
	205
	206	Int_t ntrk3d = 0;
	207	Bool_t is3d=kFALSE;
	208	if(vtx3d){
	209
	210	if(vtx3d->IsFromVertexer3D()) is3d=kTRUE;
	211	ntrk3d = vtx3d->GetNContributors();
	212	if(is3d && ntrk3d>0)good3d++;
	213	x3d = vtx3d->GetXv();
	214	xerr3d=vtx3d->GetXRes();
	215	xdiff3d = 10000.*(x3d-mcVertex[0]); // microns
	216
	217	y3d = vtx3d->GetYv();
	218	yerr3d=vtx3d->GetYRes();
	219	ydiff3d = 10000.*(y3d-mcVertex[1]); // microns
	220
	221	z3d = vtx3d->GetZv();
	222	zerr3d=vtx3d->GetZRes();
	223	zdiff3d = 10000.*(z3d-mcVertex[2]); // microns
	224
	225	}
	226
	227	xnt[0]=mcVertex[0];//x
	228	xnt[1]=mcVertex[1];//y
	229	xnt[2]=mcVertex[2];//z
	230
	231	xnt[3]=zz;
	232	xnt[4]=zdiffz;
	233	xnt[5]=zresz;
	234	xnt[6]=ntrkz;
	235
	236	xnt[7]=x3d;
	237	xnt[8]=xdiff3d;
	238	xnt[9]=xerr3d;
	239	xnt[10]=y3d;
	240	xnt[11]=ydiff3d;
	241	xnt[12]=yerr3d;
	242	xnt[13]=z3d;
	243	xnt[14]=zdiff3d;
	244	xnt[15]=zerr3d;
	245	xnt[16]=ntrk3d;
	246
	247	xnt[17]=nChGen;
	248	xnt[18]=float(ntrklets);
	249	xnt[19]=float(nrecp1);
	250	xnt[20]=float(ptype);
	251	xnt[21]=float(is3d);
	252	xnt[22]=float(eventTriggered);
	253	nt->Fill(xnt);
	254	current->cd();
	255
	256	if(optdebug){
	257	printf("Vertexer3D: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tContributors=%d \n",x3d,y3d,z3d,ntrk3d);
	258	printf("VertexerZ: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tContributors=%d \n",0.,0.,zz,ntrkz);
	259	if(isMC) printf("True Pos.: \tx=%.5f \ty=%.5f \tz(cm)=%.5f \tdN/dy=%.1f\n",mcVertex[0],mcVertex[1],mcVertex[2],nChGen);
	260	printf("Multiplicity: Tracklets=%d ClustersLay1=%d\n",ntrklets,nrecp1);
	261	}
	262
	263	}
	264	fint->cd();
	265	nt->Write();
	266	fint->Close();
	267	delete fint;
	268	if(optdebug){
	269	printf("***********************************************\n");
	270	printf("Number of Z vertices found= %d\n",goodz);
	271	printf("efficiency=%f\n",float(goodz)/float(totev));
	272	printf("Number of 3D vertices found= %d\n",good3d);
	273	printf("efficiency=%f\n",float(good3d)/float(totev));
	274	}
	275	return kTRUE;
	276	}