Interface to AliMillePede2 for ITS alignment + related classes (R. Shaoyan)
[u/mrichter/AliRoot.git] / ITS / DoVerticesSPD.C
CommitLineData
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
30Bool_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
f5826889 71 Int_t totev=runLoader->GetNumberOfEvents();
72 if(optdebug) printf("Number of events= %d\n",totev);
73
fdc48e5f 74 AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
75
f5826889 76 Double_t xnom=0.,ynom=0.;
308c2f7c 77 AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom);
78 vertz->Init("default");
79 AliITSVertexer3D *vert3d = new AliITSVertexer3D();
80 vert3d->Init("default");
792b611a 81 vert3d->SetWideFiducialRegion(40.,1.);
6b4d9537 82 vert3d->SetPileupAlgo(1);
792b611a 83 vert3d->PrintStatus();
fdc48e5f 84 vertz->SetDetTypeRec(detTypeRec);
85 vert3d->SetDetTypeRec(detTypeRec);
ecce5370 86
87 /* uncomment these lines to use diamond constrain */
88// Double_t posdiam[3]={0.03,0.1,0.};
89// Double_t sigdiam[3]={0.01,0.01,10.0};
90// AliESDVertex* diam=new AliESDVertex(posdiam,sigdiam);
91// vertz->SetVtxStart(diam);
92// vert3d->SetVtxStart(diam);
93 /* end lines to be uncommented to use diamond constrain */
f5826889 94
95 Int_t goodz=0,good3d=0;
96
97 for (Int_t iEvent = 0; iEvent < totev; iEvent++) {
98 TArrayF mcVertex(3);
99 runLoader->GetEvent(iEvent);
100 runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
101 AliGenPythiaEventHeader *evh=(AliGenPythiaEventHeader*)runLoader->GetHeader()->GenEventHeader();
102 Int_t ptype = evh->ProcessType();
103 if(optdebug){
104 printf("==============================================================\n");
105 printf("\nEvent: %d ---- Process Type = %d \n",iEvent,ptype);
106 }
107
108 AliStack* stack = runLoader->Stack();
109 TTree *treek=(TTree*)runLoader->TreeK();
110 Int_t npart = (Int_t)treek->GetEntries();
111 if(optdebug) printf("particles %d\n",npart);
112
113 Double_t dNchdy = 0.;
114
ecce5370 115 // loop on particles to get generated dN/dy
116 for(Int_t iPart=0; iPart<npart; iPart++) {
117 if(!stack->IsPhysicalPrimary(iPart)) continue;
118 TParticle* part = (TParticle*)stack->Particle(iPart);
119 if(part->GetPDG()->Charge() == 0) continue;
120 Double_t eta=part->Eta();
f5826889 121
ecce5370 122 if(TMath::Abs(eta)<1.5) dNchdy+=1.;
f5826889 123 }
124 if(optdebug) printf(" dNch/dy = %f\n",dNchdy);
125
308c2f7c 126 TTree* cltree = ITSloader->TreeR();
127
128 AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree);
ecce5370 129 AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree);
ecce5370 130 AliMultiplicity *alimult = vert3d->GetMultiplicity();
792b611a 131 Int_t ntrklets=0;
132 Int_t nrecp1=0;
133 if(alimult){
134 ntrklets=alimult->GetNumberOfTracklets() ;
135 nrecp1=ntrklets+alimult->GetNumberOfSingleClusters();
f5826889 136 }
792b611a 137
ecce5370 138
f5826889 139
140 TDirectory *current = gDirectory;
141 fint->cd();
792b611a 142 Float_t xnt[22];
f5826889 143
144 Double_t zz = 0.;
145 Double_t zresz = 0.;
146 Double_t zdiffz = 0.;
147
148 Int_t ntrkz = 0;
149 if(vtxz){
150
151 ntrkz = vtxz->GetNContributors();
152 if(ntrkz>0)goodz++;
153 zz=vtxz->GetZv();
154 zresz =vtxz->GetZRes(); // microns
155 zdiffz = 10000.*(zz-mcVertex[2]); // microns
156 }
157
158 Double_t x3d = 0.;
159 Double_t xerr3d = 0.;
160 Double_t xdiff3d = 0.;
161
162 Double_t y3d = 0.;
163 Double_t yerr3d = 0.;
164 Double_t ydiff3d = 0.;
165
166 Double_t z3d = 0.;
167 Double_t zerr3d = 0.;
168 Double_t zdiff3d = 0.;
169
170 Int_t ntrk3d = 0;
6b4d9537 171 Bool_t is3d=kFALSE;
f5826889 172 if(vtx3d){
6b4d9537 173 if(vtx3d->IsFromVertexer3D()) is3d=kTRUE;
f5826889 174 ntrk3d = vtx3d->GetNContributors();
792b611a 175 if(is3d && ntrk3d>0)good3d++;
f5826889 176 x3d = vtx3d->GetXv();
177 xerr3d=vtx3d->GetXRes();
178 xdiff3d = 10000.*(x3d-mcVertex[0]); // microns
179
180 y3d = vtx3d->GetYv();
181 yerr3d=vtx3d->GetYRes();
182 ydiff3d = 10000.*(y3d-mcVertex[1]); // microns
183
184 z3d = vtx3d->GetZv();
185 zerr3d=vtx3d->GetZRes();
186 zdiff3d = 10000.*(z3d-mcVertex[2]); // microns
187
188 }
189
190 xnt[0]=mcVertex[0];//x
191 xnt[1]=mcVertex[1];//y
192 xnt[2]=mcVertex[2];//z
193
194 xnt[3]=zz;
195 xnt[4]=zdiffz;
196 xnt[5]=zresz;
197 xnt[6]=ntrkz;
198
199 xnt[7]=x3d;
200 xnt[8]=xdiff3d;
201 xnt[9]=xerr3d;
8686a665 202 xnt[10]=y3d;
f5826889 203 xnt[11]=ydiff3d;
204 xnt[12]=yerr3d;
205 xnt[13]=z3d;
206 xnt[14]=zdiff3d;
207 xnt[15]=zerr3d;
208 xnt[16]=ntrk3d;
209
210 xnt[17]=dNchdy;
211 xnt[18]=float(ntrklets);
212 xnt[19]=float(nrecp1);
213 xnt[20]=float(ptype);
792b611a 214 xnt[21]=float(is3d);
f5826889 215 nt->Fill(xnt);
216 current->cd();
217
218 if(optdebug){
219 printf("\nVertexerZ: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",zz,ntrkz,mcVertex[2],zdiffz);
220 printf("Vertexer3D: \tz(cm)=%9.5f \tTracklets=%d \tztrue(cm)=%9.5f \tzdiff(um)=%8.2f\n",z3d,ntrk3d,mcVertex[2],zdiff3d);
221 }
222
223 }
224 fint->cd();
225 nt->Write();
226 fint->Close();
227 delete fint;
228 if(optdebug){
229 printf("***********************************************\n");
230 printf("Number of Z vertices found= %d\n",goodz);
231 printf("efficiency=%f\n",float(goodz)/float(totev));
232 printf("Number of 3D vertices found= %d\n",good3d);
233 printf("efficiency=%f\n",float(good3d)/float(totev));
234 }
235 return kTRUE;
236}