]>
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" | |
4c00c36c | 15 | #include "AliITSsegmentationSPD.h" |
f5826889 | 16 | #include "AliVertexerTracks.h" |
17 | #include "AliCDBManager.h" | |
18 | #include "AliGeomManager.h" | |
332910f8 | 19 | #include "AliGRPManager.h" |
fdc48e5f | 20 | #include "AliITSDetTypeRec.h" |
f5826889 | 21 | #include "AliITSVertexer3D.h" |
22 | #include "AliITSVertexerZ.h" | |
23 | #include "AliESDVertex.h" | |
7d766abb | 24 | #include "AliESDEvent.h" |
f5826889 | 25 | #include "AliRun.h" |
26 | #include "AliRunLoader.h" | |
27 | #include "AliGenEventHeader.h" | |
28 | #include "AliStack.h" | |
29 | #endif | |
30 | ||
31 | /* $Id$ */ | |
32 | ||
1ff24d0a | 33 | Bool_t DoVerticesSPD(Bool_t isMC=kFALSE, Int_t pileupalgo=1, Int_t optdebug=0){ |
792b611a | 34 | |
f5826889 | 35 | TFile *fint = new TFile("VertexSPD.root","recreate"); |
7d766abb | 36 | 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:isTriggered"); |
f5826889 | 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"); | |
162637e4 | 45 | man->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); |
4c00c36c | 46 | // man->SetSpecificStorage("GRP/GRP/Data",Form("local://%s",gSystem->pwd())); |
f5826889 | 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 | } | |
3714ff25 | 57 | AliGeomManager::ApplyAlignObjsFromCDB("ITS"); |
f5826889 | 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 | } | |
1ff24d0a | 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(); | |
f5826889 | 75 | } |
f5826889 | 76 | AliITSLoader* ITSloader = (AliITSLoader*) runLoader->GetLoader("ITSLoader"); |
77 | ITSloader->LoadRecPoints("read"); | |
78 | ||
3714ff25 | 79 | Int_t totev=runLoader->GetNumberOfEvents(); |
f5826889 | 80 | if(optdebug) printf("Number of events= %d\n",totev); |
81 | ||
5951029f | 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); | |
332910f8 | 95 | AliGRPManager * grpMan=new AliGRPManager(); |
96 | grpMan->ReadGRPEntry(); | |
97 | grpMan->SetMagField(); | |
98 | printf("Magnetic field set to %f T\n",AliTracker::GetBz()/10.); | |
5951029f | 99 | |
100 | AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec(); | |
4c00c36c | 101 | AliITSsegmentation* seg = new AliITSsegmentationSPD(); |
102 | detTypeRec->SetSegmentationModel(0,seg); | |
103 | ||
f5826889 | 104 | Double_t xnom=0.,ynom=0.; |
308c2f7c | 105 | AliITSVertexerZ *vertz = new AliITSVertexerZ(xnom,ynom); |
106 | vertz->Init("default"); | |
1ff24d0a | 107 | AliITSVertexer3D *vert3d = new AliITSVertexer3D(); |
308c2f7c | 108 | vert3d->Init("default"); |
792b611a | 109 | vert3d->SetWideFiducialRegion(40.,1.); |
7d766abb | 110 | vert3d->SetPileupAlgo(pileupalgo); |
792b611a | 111 | vert3d->PrintStatus(); |
fdc48e5f | 112 | vertz->SetDetTypeRec(detTypeRec); |
113 | vert3d->SetDetTypeRec(detTypeRec); | |
1ff24d0a | 114 | vert3d->SetComputeMultiplicity(kTRUE); |
ecce5370 | 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 */ | |
f5826889 | 122 | |
123 | Int_t goodz=0,good3d=0; | |
7d766abb | 124 | // Trigger mask |
125 | ULong64_t spdFO = (1 << 14); | |
126 | ULong64_t v0left = (1 << 11); | |
127 | ULong64_t v0right = (1 << 12); | |
f5826889 | 128 | |
129 | for (Int_t iEvent = 0; iEvent < totev; iEvent++) { | |
130 | TArrayF mcVertex(3); | |
131 | runLoader->GetEvent(iEvent); | |
1ff24d0a | 132 | Double_t dNchdy = 0.; |
133 | Int_t ptype = 0; | |
f5826889 | 134 | if(optdebug){ |
135 | printf("==============================================================\n"); | |
1ff24d0a | 136 | printf("Event: %d \n",iEvent); |
f5826889 | 137 | } |
1ff24d0a | 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) dNchdy+=1.; | |
156 | } | |
157 | if(optdebug) printf(" dNch/dy = %f\n",dNchdy); | |
f5826889 | 158 | } |
7d766abb | 159 | tree->GetEvent(iEvent); |
308c2f7c | 160 | |
7d766abb | 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))); | |
308c2f7c | 165 | AliESDVertex* vtxz = vertz->FindVertexForCurrentEvent(cltree); |
ecce5370 | 166 | AliESDVertex* vtx3d = vert3d->FindVertexForCurrentEvent(cltree); |
ecce5370 | 167 | AliMultiplicity *alimult = vert3d->GetMultiplicity(); |
792b611a | 168 | Int_t ntrklets=0; |
169 | Int_t nrecp1=0; | |
170 | if(alimult){ | |
171 | ntrklets=alimult->GetNumberOfTracklets() ; | |
172 | nrecp1=ntrklets+alimult->GetNumberOfSingleClusters(); | |
f5826889 | 173 | } |
ecce5370 | 174 | |
f5826889 | 175 | |
176 | TDirectory *current = gDirectory; | |
177 | fint->cd(); | |
7d766abb | 178 | Float_t xnt[23]; |
f5826889 | 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; | |
6b4d9537 | 207 | Bool_t is3d=kFALSE; |
f5826889 | 208 | if(vtx3d){ |
7d766abb | 209 | |
6b4d9537 | 210 | if(vtx3d->IsFromVertexer3D()) is3d=kTRUE; |
f5826889 | 211 | ntrk3d = vtx3d->GetNContributors(); |
792b611a | 212 | if(is3d && ntrk3d>0)good3d++; |
f5826889 | 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; | |
8686a665 | 239 | xnt[10]=y3d; |
f5826889 | 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]=dNchdy; | |
248 | xnt[18]=float(ntrklets); | |
249 | xnt[19]=float(nrecp1); | |
250 | xnt[20]=float(ptype); | |
792b611a | 251 | xnt[21]=float(is3d); |
7d766abb | 252 | xnt[22]=float(eventTriggered); |
f5826889 | 253 | nt->Fill(xnt); |
254 | current->cd(); | |
255 | ||
256 | if(optdebug){ | |
1ff24d0a | 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],dNchdy); | |
260 | printf("Multiplicity: Tracklets=%d ClustersLay1=%d\n",ntrklets,nrecp1); | |
f5826889 | 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 | } |