FXS File ID renamed to match DA output file
[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"
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
29Bool_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}