]>
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 | } |