]>
Commit | Line | Data |
---|---|---|
88ce87da | 1 | /**************************************************************************** |
2 | * This macro is supposed to do reconstruction in the ITS via Kalman * | |
3 | * tracker V2. The ITStracker is feeded with parametrized TPC tracks * | |
4 | * * | |
5 | * It does the following steps: * | |
6 | * 1) TPC tracking parameterization * | |
056891c0 | 7 | * 2) Fast points in ITS * |
8 | * 3) ITS cluster finding V2 * | |
9 | * 4) Determine z position of primary vertex * | |
10 | * - read from event header in PbPb events * | |
11 | * - determined using points on SPD in pp events (M.Masera) * | |
12 | * 5) ITS track finding V2 * | |
13 | * * in pp, redetermine the z position of the primary vertex * | |
b2bca9d4 | 14 | * using the reconstructed tracks |
056891c0 | 15 | * 6) Create a reference file with simulation info (p,PDG...) * |
88ce87da | 16 | * * |
056891c0 | 17 | * (Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it * |
18 | * from AliTPCtest.C & AliITStestV2.C by I.Belikov * | |
88ce87da | 19 | ****************************************************************************/ |
b2bca9d4 | 20 | |
056891c0 | 21 | #if !defined(__CINT__) || defined(__MAKECINT__) |
22 | //-- --- standard headers------------- | |
b2bca9d4 | 23 | #include "Riostream.h" |
056891c0 | 24 | //--------Root headers --------------- |
70521312 | 25 | #include <TFile.h> |
26 | #include <TStopwatch.h> | |
27 | #include <TObject.h> | |
056891c0 | 28 | #include <TSystem.h> |
29 | #include <TH1.h> | |
30 | #include <TH2.h> | |
31 | #include <TProfile.h> | |
32 | #include <TTree.h> | |
33 | #include <TBranch.h> | |
34 | #include <TParticle.h> | |
35 | #include <TRandom.h> | |
36 | //----- AliRoot headers --------------- | |
70521312 | 37 | #include "alles.h" |
38 | #include "AliRun.h" | |
39 | #include "AliHeader.h" | |
40 | #include "AliGenEventHeader.h" | |
41 | #include "AliMagF.h" | |
42 | #include "AliModule.h" | |
43 | #include "AliArrayI.h" | |
44 | #include "AliDigits.h" | |
45 | #include "AliITS.h" | |
46 | #include "AliTPC.h" | |
47 | #include "AliITSgeom.h" | |
48 | #include "AliITSRecPoint.h" | |
49 | #include "AliITSclusterV2.h" | |
056891c0 | 50 | #include "AliITSsimulationFastPoints.h" |
70521312 | 51 | #include "AliITStrackerV2.h" |
52 | #include "AliKalmanTrack.h" | |
53 | #include "AliTPCtrackerParam.h" | |
056891c0 | 54 | //------------------------------------- |
88ce87da | 55 | #endif |
56 | ||
056891c0 | 57 | // structure for track references |
88ce87da | 58 | typedef struct { |
59 | Int_t lab; | |
60 | Int_t pdg; | |
70521312 | 61 | Int_t mumlab; |
88ce87da | 62 | Int_t mumpdg; |
63 | Float_t Vx,Vy,Vz; | |
64 | Float_t Px,Py,Pz; | |
65 | } RECTRACK; | |
66 | ||
056891c0 | 67 | //===== Functions definition ================================================= |
68 | ||
b2bca9d4 | 69 | Int_t MarkEvtsToSkip(const Char_t *evtsName, |
70 | Bool_t *skipEvt); | |
056891c0 | 71 | |
b2bca9d4 | 72 | Int_t TPCParamTracks(const Char_t *galiceName, |
73 | const Char_t *outName, | |
74 | const Int_t coll, | |
75 | const Double_t Bfield, | |
76 | Int_t n); | |
056891c0 | 77 | |
b2bca9d4 | 78 | Int_t ITSHits2FastRecPoints(const Char_t *galiceName, |
79 | Bool_t *skipEvt, | |
80 | Int_t n); | |
056891c0 | 81 | |
b2bca9d4 | 82 | Int_t ITSFindClustersV2(const Char_t *galiceName, |
83 | const Char_t *outName, | |
84 | Bool_t *skipEvt, | |
85 | Int_t n); | |
056891c0 | 86 | |
b2bca9d4 | 87 | Int_t ZvtxFromHeader(const Char_t *galiceName, |
88 | Bool_t *skipEvt, | |
89 | Int_t n); | |
056891c0 | 90 | |
b2bca9d4 | 91 | Int_t ZvtxFromSPD(const Char_t *galiceName, |
92 | Bool_t *skipEvt, | |
93 | Int_t n); | |
056891c0 | 94 | |
b2bca9d4 | 95 | Int_t ZvtxFromTracks(const Char_t *trkame, |
96 | Bool_t *skipEvt, | |
97 | Int_t n); | |
056891c0 | 98 | |
b2bca9d4 | 99 | Int_t ZvtxFastpp(const Char_t *galiceName, |
100 | Bool_t *skipEvt, | |
101 | Int_t n); | |
056891c0 | 102 | |
b2bca9d4 | 103 | void EvalZ(TH1F *hist, |
104 | Int_t sepa, | |
105 | Float_t &av, | |
106 | Float_t &sig, | |
107 | Int_t ncoinc, | |
108 | TArrayF *zval, | |
109 | ofstream *deb); | |
88ce87da | 110 | |
b2bca9d4 | 111 | Bool_t VtxTrack(Double_t pt, |
112 | Double_t d0rphi); | |
056891c0 | 113 | |
b2bca9d4 | 114 | Double_t d0zRes(Double_t pt); |
056891c0 | 115 | |
b2bca9d4 | 116 | Int_t ITSFindTracksV2(const Char_t *galiceName, |
117 | const Char_t *inName, | |
118 | const Char_t *inName2, | |
119 | const Char_t *outName, | |
120 | Bool_t *skipEvt, | |
121 | Option_t *vtxMode, | |
122 | Int_t n); | |
056891c0 | 123 | |
b2bca9d4 | 124 | Int_t ITSMakeRefFile(const Char_t *galiceName, |
125 | const Char_t *inName, | |
126 | const Char_t *outName, | |
127 | Bool_t *skipEvt, | |
128 | Int_t n); | |
056891c0 | 129 | |
b2bca9d4 | 130 | void WriteZvtx(const Char_t *name, |
131 | Double_t *zvtx, | |
132 | Int_t n); | |
056891c0 | 133 | |
b2bca9d4 | 134 | void ReadZvtx(const Char_t *name, |
135 | Double_t *zvtx, | |
136 | Int_t n); | |
056891c0 | 137 | |
138 | //============================================================================= | |
139 | ||
b2bca9d4 | 140 | Int_t AliBarrelRec_TPCparam(Int_t n=1) { |
70521312 | 141 | |
142 | const Char_t *name=" AliBarrelRec_TPCparam"; | |
143 | cerr<<'\n'<<name<<"...\n"; | |
144 | gBenchmark->Start(name); | |
145 | ||
b2bca9d4 | 146 | const Char_t *evtsName = "EvtsToSkip.dat"; |
056891c0 | 147 | const Char_t *TPCtrkNameS = "AliTPCtracksParam.root"; |
b2bca9d4 | 148 | const Char_t *galiceName = "galice.root"; |
056891c0 | 149 | const Char_t *ITSclsName = "AliITSclustersV2.root"; |
150 | const Char_t *ITStrkName = "AliITStracksV2.root"; | |
151 | const Char_t *ITStrkName2 = "AliITStracksV2_2.root"; | |
152 | const Char_t *ITSrefName = "ITStracksRefFile.root"; | |
88ce87da | 153 | |
154 | // set here the code for the type of collision (needed for TPC tracking | |
155 | // parameterization). available collisions: | |
156 | // | |
157 | // coll = 0 -> PbPb6000 (HIJING with b<2fm) | |
056891c0 | 158 | // coll = 1 -> pp |
159 | const Int_t collcode = 1; | |
160 | const Bool_t slowVtx = kFALSE; | |
161 | const Bool_t retrack = kFALSE; | |
88ce87da | 162 | // set here the value of the magnetic field |
163 | const Double_t BfieldValue = 0.4; | |
b2bca9d4 | 164 | |
056891c0 | 165 | AliKalmanTrack::SetConvConst(100/0.299792458/BfieldValue); |
88ce87da | 166 | |
b2bca9d4 | 167 | Bool_t *skipEvt = new Bool_t[n]; |
056891c0 | 168 | |
b2bca9d4 | 169 | // Mark events that have to be skipped (read from file ascii evtsName) |
170 | for(Int_t i=0;i<n;i++) skipEvt[i] = kFALSE; | |
171 | if(!gSystem->AccessPathName(evtsName,kFileExists)) { | |
172 | MarkEvtsToSkip(evtsName,skipEvt); | |
056891c0 | 173 | } |
b2bca9d4 | 174 | |
175 | // ********** Build TPC tracks with parameterization *********** // | |
176 | TPCParamTracks(galiceName,TPCtrkNameS,collcode,BfieldValue,n); | |
177 | ||
178 | ||
179 | // ********** ITS RecPoints ************************************ // | |
180 | ITSHits2FastRecPoints(galiceName,skipEvt,n); | |
181 | ||
182 | ||
056891c0 | 183 | // ********** Find ITS clusters ******************************** // |
b2bca9d4 | 184 | ITSFindClustersV2(galiceName,ITSclsName,skipEvt,n); |
185 | ||
056891c0 | 186 | |
187 | // ********** Tracking in ITS ********************************** // | |
b2bca9d4 | 188 | Char_t *vtxMode; |
189 | // Pb-Pb | |
190 | if(collcode==0) { | |
191 | ZvtxFromHeader(galiceName,skipEvt,n); | |
192 | vtxMode="Header"; | |
193 | ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,skipEvt,vtxMode,n); | |
194 | } | |
195 | // pp | |
196 | if(collcode==1 && slowVtx) { | |
197 | ZvtxFromSPD(galiceName,skipEvt,n); | |
198 | vtxMode="SPD"; | |
199 | ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,skipEvt,vtxMode,n); | |
200 | ZvtxFromTracks(ITStrkName,skipEvt,n); | |
201 | if(retrack) { | |
202 | vtxMode="Tracks"; | |
203 | ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName2,skipEvt,vtxMode,n); | |
056891c0 | 204 | } |
b2bca9d4 | 205 | } |
206 | if(collcode==1 && !slowVtx) { | |
207 | ZvtxFastpp(galiceName,skipEvt,n); | |
208 | vtxMode="Fast"; | |
209 | ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,skipEvt,vtxMode,n); | |
056891c0 | 210 | } |
211 | ||
056891c0 | 212 | // ********** Make ITS tracks reference file ******************* // |
213 | ITSMakeRefFile(galiceName,ITStrkName,ITSrefName,skipEvt,n); | |
b2bca9d4 | 214 | |
215 | ||
056891c0 | 216 | |
217 | delete [] skipEvt; | |
218 | ||
70521312 | 219 | gBenchmark->Stop(name); |
220 | gBenchmark->Show(name); | |
88ce87da | 221 | |
b2bca9d4 | 222 | return 0; |
056891c0 | 223 | } |
224 | //----------------------------------------------------------------------------- | |
b2bca9d4 | 225 | Int_t MarkEvtsToSkip(const Char_t *evtsName,Bool_t *skipEvt) { |
88ce87da | 226 | |
056891c0 | 227 | cerr<<"\n*******************************************************************\n"; |
228 | cerr<<"\nChecking for events to skip...\n"; | |
229 | ||
230 | Int_t evt,ncol; | |
231 | ||
232 | FILE *f = fopen(evtsName,"r"); | |
233 | while(1) { | |
234 | ncol = fscanf(f,"%d",&evt); | |
235 | if(ncol<1) break; | |
b2bca9d4 | 236 | skipEvt[evt] = kTRUE; |
056891c0 | 237 | cerr<<" ev. "<<evt<<" will be skipped\n"; |
238 | } | |
239 | fclose(f); | |
88ce87da | 240 | |
056891c0 | 241 | return 0; |
242 | } | |
243 | //----------------------------------------------------------------------------- | |
244 | Int_t TPCParamTracks(const Char_t *galiceName,const Char_t *outName, | |
88ce87da | 245 | const Int_t coll,const Double_t Bfield,Int_t n) { |
246 | ||
247 | cerr<<"\n*******************************************************************\n"; | |
88ce87da | 248 | |
249 | const Char_t *name="TPCParamTracks"; | |
250 | cerr<<'\n'<<name<<"...\n"; | |
251 | gBenchmark->Start(name); | |
252 | ||
056891c0 | 253 | TFile *outFile=TFile::Open(outName,"recreate"); |
254 | TFile *inFile =TFile::Open(galiceName); | |
255 | ||
256 | AliTPCtrackerParam tracker(coll,Bfield,n); | |
257 | tracker.BuildTPCtracks(inFile,outFile); | |
88ce87da | 258 | |
259 | delete gAlice; gAlice=0; | |
260 | ||
056891c0 | 261 | inFile->Close(); |
262 | outFile->Close(); | |
88ce87da | 263 | |
264 | gBenchmark->Stop(name); | |
265 | gBenchmark->Show(name); | |
266 | ||
056891c0 | 267 | return 0; |
88ce87da | 268 | } |
056891c0 | 269 | //----------------------------------------------------------------------------- |
b2bca9d4 | 270 | Int_t ITSHits2FastRecPoints(const Char_t *galiceName,Bool_t *skipEvt,Int_t n) { |
056891c0 | 271 | |
272 | cerr<<"\n*******************************************************************\n"; | |
273 | ||
274 | const Char_t *name="ITSHits2FastRecPoints"; | |
275 | cerr<<'\n'<<name<<"...\n"; | |
276 | gBenchmark->Start(name); | |
277 | ||
278 | Int_t nsignal=25; | |
279 | Int_t size=-1; | |
88ce87da | 280 | |
056891c0 | 281 | // Connect the Root Galice file containing Geometry, Kine and Hits |
282 | TFile *file = TFile::Open(galiceName,"UPDATE"); | |
88ce87da | 283 | |
056891c0 | 284 | gAlice = (AliRun*)file->Get("gAlice"); |
285 | ||
286 | ||
287 | AliITS *ITS = (AliITS*)gAlice->GetModule("ITS"); | |
288 | if(!ITS) return 1; | |
289 | ||
290 | // Set the simulation model | |
291 | for(Int_t i=0;i<3;i++) { | |
292 | ITS->SetSimulationModel(i,new AliITSsimulationFastPoints()); | |
293 | } | |
294 | ||
295 | ||
296 | // | |
297 | // Event Loop | |
298 | // | |
299 | for(Int_t ev=0; ev<n; ev++) { | |
b2bca9d4 | 300 | if(skipEvt[ev]) continue; |
056891c0 | 301 | cerr<<" --- Processing event "<<ev<<" ---"<<endl; |
302 | Int_t nparticles = gAlice->GetEvent(ev); | |
303 | cerr<<"Number of particles: "<<nparticles<<endl; | |
304 | gAlice->SetEvent(ev); | |
305 | if(!gAlice->TreeR()) gAlice->MakeTree("R"); | |
306 | if(!gAlice->TreeR()->GetBranch("ITSRecPointsF")) { | |
307 | ITS->MakeBranch("RF"); | |
308 | if(nparticles <= 0) return 1; | |
309 | ||
310 | Int_t bgr_ev=Int_t(ev/nsignal); | |
311 | //printf("bgr_ev %d\n",bgr_ev); | |
312 | ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," "); | |
313 | } | |
314 | } // event loop | |
315 | ||
316 | delete gAlice; gAlice=0; | |
317 | file->Close(); | |
318 | ||
319 | gBenchmark->Stop(name); | |
320 | gBenchmark->Show(name); | |
321 | ||
322 | return 0; | |
323 | } | |
324 | //----------------------------------------------------------------------------- | |
325 | Int_t ITSFindClustersV2(const Char_t *galiceName,const Char_t *outName, | |
b2bca9d4 | 326 | Bool_t *skipEvt,Int_t n) { |
056891c0 | 327 | // |
328 | // This function converts AliITSRecPoint(s) to AliITSclusterV2 | |
329 | // | |
88ce87da | 330 | cerr<<"\n*******************************************************************\n"; |
331 | ||
056891c0 | 332 | const Char_t *name="ITSFindClustersV2"; |
88ce87da | 333 | cerr<<'\n'<<name<<"...\n"; |
334 | gBenchmark->Start(name); | |
70521312 | 335 | |
70521312 | 336 | |
056891c0 | 337 | TFile *inFile=TFile::Open(galiceName); |
338 | if(!inFile->IsOpen()) { cerr<<"Can't open galice.root !\n"; return 1; } | |
339 | ||
340 | if(!(gAlice=(AliRun*)inFile->Get("gAlice"))) { | |
341 | cerr<<"Can't find gAlice !\n"; | |
88ce87da | 342 | return 1; |
343 | } | |
88ce87da | 344 | AliITS *ITS = (AliITS*)gAlice->GetModule("ITS"); |
056891c0 | 345 | if(!ITS) { cerr<<"Can't find the ITS !\n"; return 1; } |
88ce87da | 346 | AliITSgeom *geom=ITS->GetITSgeom(); |
056891c0 | 347 | |
348 | TFile *outFile = TFile::Open(outName,"recreate"); | |
88ce87da | 349 | geom->Write(); |
88ce87da | 350 | |
056891c0 | 351 | // loop on events |
352 | for(Int_t ev=0; ev<n; ev++) { | |
056891c0 | 353 | if(skipEvt[ev]) continue; |
056891c0 | 354 | cerr<<"--- Processing event "<<ev<<" ---"<<endl; |
355 | inFile->cd(); | |
356 | gAlice->GetEvent(ev); | |
357 | ||
358 | TClonesArray *clusters = new TClonesArray("AliITSclusterV2",10000); | |
359 | Char_t cname[100]; | |
88ce87da | 360 | sprintf(cname,"TreeC_ITS_%d",ev); |
056891c0 | 361 | TTree *cTree = new TTree(cname,"ITS clusters"); |
88ce87da | 362 | cTree->Branch("Clusters",&clusters); |
056891c0 | 363 | |
364 | TTree *pTree=gAlice->TreeR(); | |
365 | if(!pTree) { cerr<<"Can't get TreeR !\n"; return 1; } | |
366 | TBranch *branch = pTree->GetBranch("ITSRecPointsF"); | |
367 | if(!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1; } | |
368 | TClonesArray *points = new TClonesArray("AliITSRecPoint",10000); | |
88ce87da | 369 | branch->SetAddress(&points); |
056891c0 | 370 | |
88ce87da | 371 | TClonesArray &cl=*clusters; |
372 | Int_t nclusters=0; | |
056891c0 | 373 | Int_t nentr=(Int_t)branch->GetEntries(); |
88ce87da | 374 | |
056891c0 | 375 | cerr<<"Number of entries in TreeR_RF: "<<nentr<<endl; |
376 | ||
b2bca9d4 | 377 | Float_t *lp = new Float_t[5]; |
378 | Int_t *lab = new Int_t[6]; | |
056891c0 | 379 | |
380 | for(Int_t i=0; i<nentr; i++) { | |
381 | points->Clear(); | |
382 | branch->GetEvent(i); | |
383 | Int_t ncl=points->GetEntriesFast(); if(ncl==0){cTree->Fill();continue;} | |
88ce87da | 384 | Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det); |
385 | Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); | |
386 | Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot); | |
387 | Double_t yshift = x*rot[0] + y*rot[1]; | |
388 | Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1); | |
88ce87da | 389 | nclusters+=ncl; |
056891c0 | 390 | |
391 | Float_t kmip=1; // ADC->mip normalization factor for the SDD and SSD | |
b2bca9d4 | 392 | if(lay==4 || lay==3){kmip=280.;}; |
393 | if(lay==6 || lay==5){kmip=38.;}; | |
056891c0 | 394 | |
395 | for(Int_t j=0; j<ncl; j++) { | |
88ce87da | 396 | AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j); |
056891c0 | 397 | lp[0]=-p->GetX()-yshift; if(lay==1) lp[0]=-lp[0]; |
88ce87da | 398 | lp[1]=p->GetZ()+zshift; |
399 | lp[2]=p->GetSigmaX2(); | |
400 | lp[3]=p->GetSigmaZ2(); | |
056891c0 | 401 | lp[4]=p->GetQ(); lp[4]/=kmip; |
88ce87da | 402 | lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2); |
403 | lab[3]=ndet; | |
b2bca9d4 | 404 | |
88ce87da | 405 | Int_t label=lab[0]; |
056891c0 | 406 | if(label>=0) { |
407 | TParticle *part=(TParticle*)gAlice->Particle(label); | |
408 | label=-3; | |
409 | while (part->P() < 0.005) { | |
410 | Int_t m=part->GetFirstMother(); | |
411 | if(m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;} | |
412 | label=m; | |
413 | part=(TParticle*)gAlice->Particle(label); | |
414 | } | |
415 | if (lab[1]<0) lab[1]=label; | |
416 | else if(lab[2]<0) lab[2]=label; | |
417 | else cerr<<"No empty labels !\n"; | |
88ce87da | 418 | } |
b2bca9d4 | 419 | |
88ce87da | 420 | new(cl[j]) AliITSclusterV2(lab,lp); |
421 | } | |
422 | cTree->Fill(); clusters->Delete(); | |
423 | points->Delete(); | |
424 | } | |
056891c0 | 425 | |
426 | outFile->cd(); | |
88ce87da | 427 | cTree->Write(); |
056891c0 | 428 | |
88ce87da | 429 | cerr<<"Number of clusters: "<<nclusters<<endl; |
b2bca9d4 | 430 | delete [] lp; |
431 | delete [] lab; | |
88ce87da | 432 | delete cTree; delete clusters; delete points; |
056891c0 | 433 | |
434 | } // end loop on events | |
88ce87da | 435 | |
88ce87da | 436 | |
88ce87da | 437 | delete gAlice; gAlice=0; |
056891c0 | 438 | |
439 | inFile->Close(); | |
440 | outFile->Close(); | |
441 | ||
88ce87da | 442 | gBenchmark->Stop(name); |
443 | gBenchmark->Show(name); | |
88ce87da | 444 | |
056891c0 | 445 | return 0; |
446 | } | |
447 | //----------------------------------------------------------------------------- | |
b2bca9d4 | 448 | Int_t ZvtxFromHeader(const Char_t *galiceName, |
449 | Bool_t *skipEvt,Int_t n) { | |
88ce87da | 450 | |
88ce87da | 451 | cerr<<"\n*******************************************************************\n"; |
452 | ||
056891c0 | 453 | const Char_t *name="ZvtxFromHeader"; |
88ce87da | 454 | cerr<<'\n'<<name<<"...\n"; |
455 | gBenchmark->Start(name); | |
70521312 | 456 | |
056891c0 | 457 | TFile *galice = TFile::Open(galiceName); |
88ce87da | 458 | |
056891c0 | 459 | Double_t *zvtx = new Double_t[n]; |
460 | ||
461 | for(Int_t ev=0; ev<n; ev++){ | |
b2bca9d4 | 462 | if(skipEvt[ev]) continue; |
056891c0 | 463 | cerr<<" --- Processing event "<<ev<<" ---"<<endl; |
88ce87da | 464 | |
056891c0 | 465 | TArrayF o = 0; |
88ce87da | 466 | o.Set(3); |
88ce87da | 467 | AliHeader* header = 0; |
468 | TTree* treeE = (TTree*)gDirectory->Get("TE"); | |
469 | treeE->SetBranchAddress("Header",&header); | |
470 | treeE->GetEntry(ev); | |
471 | AliGenEventHeader* genHeader = header->GenEventHeader(); | |
472 | if(genHeader) { | |
473 | // get primary vertex position | |
474 | genHeader->PrimaryVertex(o); | |
056891c0 | 475 | // set position of primary vertex |
476 | zvtx[ev] = (Double_t)o[2]; | |
477 | } else { | |
478 | cerr<<" ! event header not found : setting z vertex to 0 !"<<endl; | |
479 | zvtx[ev] = 0.; | |
480 | } | |
481 | delete header; | |
482 | } | |
483 | ||
484 | galice->Close(); | |
485 | ||
486 | // Write vertices to file | |
b2bca9d4 | 487 | WriteZvtx("zvtxHeader.dat",zvtx,n); |
056891c0 | 488 | |
489 | delete [] zvtx; | |
490 | ||
491 | gBenchmark->Stop(name); | |
492 | gBenchmark->Show(name); | |
493 | ||
494 | return 0; | |
495 | } | |
496 | //----------------------------------------------------------------------------- | |
b2bca9d4 | 497 | Int_t ZvtxFastpp(const Char_t *galiceName, |
498 | Bool_t *skipEvt,Int_t n) { | |
056891c0 | 499 | |
500 | cerr<<"\n*******************************************************************\n"; | |
501 | ||
502 | const Char_t *name="ZvtxFastpp"; | |
503 | cerr<<'\n'<<name<<"...\n"; | |
504 | gBenchmark->Start(name); | |
505 | ||
506 | TFile *galice = TFile::Open(galiceName); | |
507 | AliRun *gAlice = (AliRun*)galice->Get("gAlice"); | |
508 | ||
509 | Double_t *zvtx = new Double_t[n]; | |
510 | ||
511 | TParticle *part; | |
512 | Int_t nPart,b; | |
513 | Double_t dNchdy,E,theta,eta,zvtxTrue,sigmaz,probVtx; | |
514 | ||
515 | for(Int_t ev=0; ev<n; ev++) { | |
b2bca9d4 | 516 | if(skipEvt[ev]) continue; |
056891c0 | 517 | cerr<<" --- Processing event "<<ev<<" ---"<<endl; |
518 | ||
519 | TArrayF o = 0; | |
520 | o.Set(3); | |
521 | AliHeader* header = 0; | |
522 | TTree* treeE = (TTree*)galice->Get("TE"); | |
523 | treeE->SetBranchAddress("Header",&header); | |
524 | treeE->GetEntry(ev); | |
525 | AliGenEventHeader* genHeader = header->GenEventHeader(); | |
526 | if(genHeader) { | |
527 | // get primary vertex position | |
528 | genHeader->PrimaryVertex(o); | |
529 | // set position of primary vertex | |
530 | zvtxTrue = (Double_t)o[2]; | |
531 | } else { | |
532 | cerr<<" ! event header not found : setting z vertex to 0 !"<<endl; | |
533 | zvtxTrue = 0.; | |
88ce87da | 534 | } |
056891c0 | 535 | delete header; |
536 | ||
537 | ||
538 | // calculate dNch/dy pf the event (charged primaries in |eta|<0.5) | |
539 | nPart = (Int_t)gAlice->GetEvent(ev); | |
540 | dNchdy = 0.; | |
541 | for(b=0; b<nPart; b++) { | |
542 | part = (TParticle*)gAlice->Particle(b); | |
543 | if(TMath::Abs(part->GetPdgCode())<10) continue; // reject partons | |
544 | if(TMath::Abs(part->Vx()*part->Vx()+part->Vy()*part->Vy())>0.01) continue; // reject secondaries | |
545 | E = part->Energy(); | |
546 | if(E>6900.) continue; // reject incoming protons | |
547 | theta = part->Theta(); | |
548 | if(theta<.1 || theta>3.) continue; | |
549 | eta = -TMath::Log(TMath::Tan(theta/2.)); | |
550 | if(TMath::Abs(eta)>0.5) continue; // count particles in |eta|<0.5 | |
551 | dNchdy+=1.; | |
552 | } | |
553 | ||
554 | // get sigma(z) corresponding to the mult of this event | |
555 | if(dNchdy>1.) { | |
556 | sigmaz = 0.0417/TMath::Sqrt(dNchdy); | |
557 | } else { | |
558 | sigmaz = 0.0500; | |
559 | } | |
560 | // smear the original position of the primary vertex | |
561 | zvtx[ev] = gRandom->Gaus(zvtxTrue,sigmaz); | |
562 | ||
563 | // compute the probability that the vertex is found | |
564 | probVtx = 1.; | |
565 | if(dNchdy<24.) probVtx = 0.85+0.00673*dNchdy; | |
566 | if(dNchdy<14.) probVtx = 0.85; | |
567 | ||
568 | if(gRandom->Rndm()>probVtx) zvtx[ev] = -1000.;// no zvtx for this event | |
569 | ||
570 | } | |
571 | galice->Close(); | |
572 | ||
573 | // Write vertices to file | |
b2bca9d4 | 574 | WriteZvtx("zvtxFastpp.dat",zvtx,n); |
056891c0 | 575 | |
576 | delete [] zvtx; | |
577 | //delete gAlice; gAlice=0; | |
578 | ||
579 | ||
580 | gBenchmark->Stop(name); | |
581 | gBenchmark->Show(name); | |
582 | ||
583 | return 0; | |
584 | } | |
585 | //----------------------------------------------------------------------------- | |
b2bca9d4 | 586 | Int_t ZvtxFromSPD(const Char_t *galiceName,Bool_t *skipEvt,Int_t n) { |
056891c0 | 587 | |
588 | ||
589 | Double_t *zvtx = new Double_t[n]; | |
590 | ||
591 | Bool_t debug = kFALSE; | |
592 | ofstream fildeb; | |
593 | if(debug)fildeb.open("FindZV.out"); | |
594 | ofstream filou; | |
595 | //filou.open("zcoor.dat"); | |
596 | const Float_t kPi=TMath::Pi(); | |
597 | const Int_t kFirstL1=0; | |
598 | const Int_t kLastL1=79; | |
599 | const Int_t kFirstL2=80; | |
600 | const Int_t kLastL2=239; | |
601 | // const Float_t kDiffPhiMax=0.0005; | |
602 | const Float_t kDiffPhiMax=0.01; | |
603 | const Float_t kphl=0.; | |
604 | const Float_t kphh=kPi/18.; | |
605 | TDirectory *currdir = gDirectory; | |
606 | TH2F *hz2z1 = 0; | |
607 | TH1F *zvdis = 0; | |
608 | TH2F *dpvsz = 0; | |
609 | TProfile *scoc = new TProfile("scoc","Zgen - Zmeas vs occ lay 1",5,10.,60.); | |
610 | TProfile *prof = new TProfile("prof","z2 vs z1",50,-15.,15.); | |
611 | TH1F *phdiff = new TH1F("phdiff","#Delta #phi",100,kphl,kphh); | |
612 | TH1F *phdifsame = new TH1F("phdifsame","#Delta #phi same track",100,kphl,kphh); | |
613 | if(gAlice){delete gAlice; gAlice = 0;} | |
614 | TFile *file = TFile::Open(galiceName); | |
615 | gAlice = (AliRun*)file->Get("gAlice"); | |
616 | ||
617 | AliITS *ITS = (AliITS*)gAlice->GetModule("ITS"); | |
618 | AliITSgeom *geom = ITS->GetITSgeom(); | |
619 | TTree *TR=0; | |
620 | Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.; | |
621 | Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.; | |
622 | Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.; | |
623 | Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.; | |
624 | char name[30]; | |
625 | Float_t avesca=0; | |
626 | Float_t avesca2=0; | |
627 | Int_t N=0; | |
628 | for(Int_t event=0; event<n; event++){ | |
b2bca9d4 | 629 | if(skipEvt[event]) continue; |
056891c0 | 630 | sprintf(name,"event_%d",event); |
631 | hz2z1=new TH2F(name,"z2 vs z1",50,-15.,15.,50,-15.,15.); | |
632 | hz2z1->SetDirectory(currdir); | |
633 | gAlice->GetEvent(event); | |
634 | TParticle * part = gAlice->Particle(0); | |
635 | Float_t oriz=part->Vz(); | |
636 | cout<<"\n ==============================================================\n"; | |
637 | cout<<" Processing event "<<event<<" "<<oriz<<endl; | |
638 | cout<<" ==============================================================\n"; | |
639 | if(debug){ | |
640 | fildeb<<"\n ==============================================================\n"; | |
641 | fildeb<<" Processing event "<<event<<" "<<oriz<<endl; | |
642 | fildeb<<" ==============================================================\n"; | |
643 | } | |
644 | file->cd(); | |
645 | TR = gAlice->TreeR(); | |
646 | if(!TR) cerr<<"TreeR not found\n"; | |
647 | TBranch *branch=TR->GetBranch("ITSRecPointsF"); | |
648 | if(!branch) cerr<<"Branch ITSRecPointsF not found\n"; | |
649 | TClonesArray *points = new TClonesArray("AliITSRecPoint",10000); | |
650 | branch->SetAddress(&points); | |
651 | Int_t nentr=(Int_t)branch->GetEntries(); | |
652 | Float_t zave=0; | |
653 | Float_t rmszav=0; | |
654 | Float_t zave2=0; | |
655 | Int_t firipixe=0; | |
656 | cout<<endl; | |
657 | for(Int_t module= kFirstL1; module<=kLastL1;module++){ | |
658 | points->Clear(); | |
659 | branch->GetEvent(module); | |
660 | Int_t nrecp1 = points->GetEntriesFast(); | |
661 | //cerr<<nrecp1<<endl; | |
662 | for(Int_t i=0; i<nrecp1;i++){ | |
663 | AliITSRecPoint *current = (AliITSRecPoint*)points->UncheckedAt(i); | |
664 | lc[0]=current->GetX(); | |
665 | lc[2]=current->GetZ(); | |
666 | geom->LtoG(module,lc,gc); | |
667 | zave+=gc[2]; | |
668 | zave2+=gc[2]*gc[2]; | |
669 | firipixe++; | |
670 | } | |
671 | points->Delete(); | |
672 | } | |
673 | if(firipixe>1){ | |
674 | rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1)); | |
675 | zave=zave/firipixe; | |
676 | cout<<"Z aver first layer = "<<zave<<" RMS= "<<rmszav<<" pixels = "<<firipixe<<endl; | |
677 | if(debug)fildeb<<"Z aver first layer = "<<zave<<" RMS= "<<rmszav<<" pixels = "<<firipixe<<endl; | |
678 | } | |
679 | else { | |
680 | cout<<"No rec points on first layer for this event\n"; | |
681 | if(debug)fildeb<<"No rec points on first layer for this event\n"; | |
682 | } | |
683 | Float_t zlim1=zave-rmszav; | |
684 | Float_t zlim2=zave+rmszav; | |
685 | Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.); | |
686 | zlim2=zlim1 + sepa/10.; | |
687 | sprintf(name,"pz_ev_%d",event); | |
688 | dpvsz=new TH2F(name,"#Delta #phi vs Z",sepa,zlim1,zlim2,100,0.,0.03); | |
689 | dpvsz->SetDirectory(currdir); | |
690 | sprintf(name,"z_ev_%d",event); | |
691 | zvdis=new TH1F(name,"zv distr",sepa,zlim1,zlim2); | |
692 | zvdis->SetDirectory(currdir); | |
693 | cout<<"Z limits: "<<zlim1<<" "<<zlim2<<endl; | |
694 | if(debug)fildeb<<"Z limits: "<<zlim1<<" "<<zlim2<<endl; | |
695 | Int_t sizarr=100; | |
696 | TArrayF *zval = new TArrayF(sizarr); | |
697 | Int_t ncoinc=0; | |
698 | for(Int_t module= kFirstL1; module<=kLastL1;module++){ | |
699 | //cout<<"processing module "<<module<<" \r"; | |
700 | branch->GetEvent(module); | |
701 | Int_t nrecp1 = points->GetEntriesFast(); | |
702 | TObjArray *poiL1 = new TObjArray(nrecp1); | |
703 | for(Int_t i=0; i<nrecp1;i++)poiL1->AddAt(points->UncheckedAt(i),i); | |
704 | //ITS->ResetRecPoints(); | |
705 | points->Delete(); | |
706 | for(Int_t i=0; i<nrecp1;i++){ | |
707 | AliITSRecPoint *current = (AliITSRecPoint*)poiL1->UncheckedAt(i); | |
708 | lc[0]=current->GetX(); | |
709 | lc[2]=current->GetZ(); | |
710 | geom->LtoG(module,lc,gc); | |
711 | Float_t r1=TMath::Sqrt(gc[0]*gc[0]+gc[1]*gc[1]); | |
712 | Float_t phi1 = TMath::ATan2(gc[1],gc[0]); | |
713 | if(phi1<0)phi1=2*kPi+phi1; | |
714 | Int_t lab1 = current->GetLabel(0); | |
715 | Float_t phi1d=phi1*180./kPi; | |
716 | //cerr<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1d<<" "<<lab1<<" \n"; | |
717 | for(Int_t modul2=kFirstL2; modul2<=kLastL2; modul2++){ | |
718 | branch->GetEvent(modul2); | |
719 | Int_t nrecp2 = points->GetEntriesFast(); | |
720 | for(Int_t j=0; j<nrecp2;j++){ | |
721 | AliITSRecPoint *recp = (AliITSRecPoint*)points->UncheckedAt(j); | |
722 | lc2[0]=recp->GetX(); | |
723 | lc2[2]=recp->GetZ(); | |
724 | geom->LtoG(modul2,lc2,gc2); | |
725 | Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]); | |
726 | Float_t zr0=(r2*gc[2]-r1*gc2[2])/(r2-r1); | |
727 | Int_t lab2 = recp->GetLabel(0); | |
728 | Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]); | |
729 | if(phi2<0)phi2=2.*kPi+phi2; | |
730 | Float_t diff = TMath::Abs(phi2-phi1); | |
731 | if(diff>kPi)diff=2.*kPi-diff; | |
732 | if(lab1==lab2)phdifsame->Fill(diff); | |
733 | if(zr0>zlim1 && zr0<zlim2){ | |
734 | phdiff->Fill(diff); | |
735 | dpvsz->Fill(zr0,diff); | |
736 | if(diff<kDiffPhiMax ){ | |
737 | zvdis->Fill(zr0); | |
738 | zval->AddAt(zr0,ncoinc); | |
739 | ncoinc++; | |
740 | if(ncoinc==(sizarr-1)){ | |
741 | sizarr+=100; | |
742 | zval->Set(sizarr); | |
743 | } | |
744 | Float_t phi2d=phi2*180./kPi; | |
745 | Float_t diffd=diff*180./kPi; | |
746 | //cerr<<"module2 "<<modul2<<" "<<gc2[0]<<" "<<gc2[1]<<" "<<gc2[2]<<" "<<phi2d<<" "<<diffd<<" "<<lab2<<" \n"; | |
747 | // cout<<"zr0= "<<zr0<<endl; | |
748 | hz2z1->Fill(gc[2],gc2[2]); | |
749 | prof->Fill(gc[2]-oriz,gc2[2]-oriz); | |
750 | } | |
751 | } | |
752 | } | |
753 | points->Delete(); | |
754 | } | |
755 | } | |
756 | delete poiL1; | |
757 | } // loop on modules | |
758 | cout<<endl; | |
759 | Float_t zcmp=0; | |
760 | Float_t zsig=0; | |
761 | EvalZ(zvdis,sepa,zcmp,zsig,ncoinc,zval,&fildeb); | |
762 | if(zcmp!=-100 && zsig!=-100){ | |
763 | N++; | |
764 | Float_t scarto = (oriz-zcmp)*10000.; | |
765 | Float_t ascarto=TMath::Abs(scarto); | |
766 | scoc->Fill(firipixe,ascarto); | |
767 | avesca+=ascarto; | |
768 | avesca2+=scarto*scarto; | |
769 | if(debug)fildeb<<"Mean value: "<<zcmp<<" +/-"<<zsig<<" Zgen- Zmeas (micron)= "<<scarto<<endl; | |
770 | cout<<"Mean value: "<<zcmp<<" RMS: "<<zsig<<" Zgen- Zmeas (micron)= "<<scarto<<endl; | |
771 | //filou<<event<<" "<<zcmp<<" "<<zsig<<" "<<oriz<<" "<<scarto<<endl; | |
772 | zvtx[event] = (Double_t)zcmp; | |
773 | } | |
774 | else { | |
775 | cout<<"Not enough points in event "<<event<<endl; | |
776 | if(debug)fildeb<<"Not enough points\n"; | |
777 | //filou<<event<<" "<<-1000.<<" "<<0.<<" "<<oriz<<" "<<0.<<endl; | |
778 | zvtx[event] = -1000.; | |
779 | } | |
780 | delete zval; | |
781 | } // loop on events | |
782 | file->Close(); | |
783 | //hz2z1->Draw("box"); | |
784 | ||
785 | // Write vertices to file | |
b2bca9d4 | 786 | WriteZvtx("zvtxSPD.dat",zvtx,n); |
056891c0 | 787 | delete [] zvtx; |
788 | ||
789 | if(N>1) { | |
790 | avesca2=TMath::Sqrt((avesca2/(N-1)-avesca*avesca/N/(N-1))/N); | |
791 | avesca=avesca/N; | |
792 | } else { | |
793 | avesca2 = 0; | |
794 | } | |
795 | cout<<"\n \n ==========================================================\n"; | |
796 | cout<<"Number of events with vertex estimate "<<N<<endl; | |
797 | cout<<"Average of the (abs) Zgen-Zmeas : "<<avesca; | |
798 | cout<<" +/- "<<avesca2<<" micron"<<endl; | |
799 | if(debug){ | |
800 | fildeb<<"\n \n ==========================================================\n"; | |
801 | fildeb<<"Number of events with vertex estimate "<<N<<endl; | |
802 | fildeb<<"Average of the (abs) Zgen-Zmeas : "<<avesca; | |
803 | fildeb<<" +/- "<<avesca2<<endl; | |
804 | } | |
805 | return 0; | |
806 | } | |
807 | ||
808 | ||
809 | //----------------------------------------------------------------------------- | |
810 | void EvalZ(TH1F *hist,Int_t sepa,Float_t &av, Float_t &sig,Int_t ncoinc, | |
811 | TArrayF *zval,ofstream *deb) { | |
812 | ||
813 | av=0; | |
814 | sig=0; | |
815 | Int_t N=0; | |
816 | Int_t NbinNotZero=0; | |
817 | for(Int_t i=1;i<=sepa;i++){ | |
818 | Float_t cont=hist->GetBinContent(i); | |
819 | av+=cont; | |
820 | sig+=cont*cont; | |
821 | N++; | |
822 | if(cont!=0)NbinNotZero++; | |
823 | } | |
824 | if(NbinNotZero==0){av=-100; sig=-100; return;} | |
825 | if(NbinNotZero==1){sig=-100; return;} | |
826 | sig=TMath::Sqrt(sig/(N-1)-av*av/N/(N-1)); | |
827 | av=av/N; | |
828 | Float_t val1=hist->GetBinLowEdge(sepa); | |
829 | Float_t val2=hist->GetBinLowEdge(1); | |
830 | for(Int_t i=1;i<=sepa;i++){ | |
831 | Float_t cont=hist->GetBinContent(i); | |
832 | if(cont>(av+sig*2.)){ | |
833 | Float_t curz=hist->GetBinLowEdge(i); | |
834 | if(curz<val1)val1=curz; | |
835 | if(curz>val2)val2=curz; | |
836 | } | |
837 | } | |
838 | val2+=hist->GetBinWidth(1); | |
839 | cout<<"Values for Z finding: "<<val1<<" "<<val2<<endl; | |
840 | if(deb)(*deb)<<"Values for Z finding: "<<val1<<" "<<val2<<endl; | |
841 | av=0; | |
842 | sig=0; | |
843 | N=0; | |
844 | for(Int_t i=0; i<ncoinc; i++){ | |
845 | Float_t z=zval->At(i); | |
846 | if(z<val1)continue; | |
847 | if(z>val2)continue; | |
848 | av+=z; | |
849 | sig+=z*z; | |
850 | N++; | |
851 | } | |
852 | if(N<=1){av=-100; sig=-100; return;} | |
853 | sig=TMath::Sqrt((sig/(N-1)-av*av/N/(N-1))/N); | |
854 | av=av/N; | |
855 | return; | |
856 | } | |
857 | //----------------------------------------------------------------------------- | |
b2bca9d4 | 858 | Int_t ZvtxFromTracks(const Char_t *trkName,Bool_t *skipEvt,Int_t n) { |
056891c0 | 859 | |
860 | cerr<<"\n*******************************************************************\n"; | |
861 | ||
862 | const Char_t *name="ZvtxFromTracks"; | |
863 | cerr<<'\n'<<name<<"...\n"; | |
864 | gBenchmark->Start(name); | |
865 | ||
866 | Double_t *zvtx = new Double_t[n]; | |
867 | Double_t *zvtxSPD = new Double_t[n]; | |
b2bca9d4 | 868 | ReadZvtx("zvtxSPD.dat",zvtxSPD,n); |
056891c0 | 869 | |
870 | TFile *itstrk = TFile::Open(trkName); | |
871 | ||
872 | ||
873 | Int_t nTrks,i,nUsedTrks; | |
874 | Double_t pt,d0rphi,d0z; | |
875 | Double_t zvtxTrks,ezvtxTrks,sumWeights; | |
876 | ||
877 | for(Int_t ev=0; ev<n; ev++){ | |
878 | cerr<<" --- Processing event "<<ev<<" ---"<<endl; | |
879 | ||
880 | // Tree with ITS tracks | |
881 | Char_t tname[100]; | |
882 | sprintf(tname,"TreeT_ITS_%d",ev); | |
883 | ||
884 | TTree *tracktree=(TTree*)itstrk->Get(tname); | |
885 | if(!tracktree) continue; | |
886 | AliITStrackV2 *itstrack=new AliITStrackV2; | |
887 | tracktree->SetBranchAddress("tracks",&itstrack); | |
888 | nTrks = (Int_t)tracktree->GetEntries(); | |
889 | nUsedTrks = 0; | |
890 | zvtxTrks = 0.; | |
891 | ezvtxTrks = 0.; | |
892 | sumWeights = 0.; | |
893 | ||
894 | //cerr<<" ITS tracks: "<<nTrks<<endl; | |
895 | // loop on tracks | |
896 | for(i=0; i<nTrks; i++) { | |
897 | tracktree->GetEvent(i); | |
898 | // get pt and impact parameters | |
899 | pt = 1./TMath::Abs(itstrack->Get1Pt()); | |
900 | d0rphi = TMath::Abs(itstrack->GetD()); | |
901 | itstrack->PropagateToVertex(); | |
902 | d0z = itstrack->GetZ(); | |
903 | // check if the track is from (0,0) in (x,y) | |
904 | if(!VtxTrack(pt,d0rphi)) continue; | |
905 | nUsedTrks++; | |
906 | zvtxTrks += d0z/d0zRes(pt)/d0zRes(pt); | |
907 | sumWeights += 1./d0zRes(pt)/d0zRes(pt); | |
908 | ||
909 | } // loop on tracks | |
910 | ||
911 | if(nUsedTrks>1) { | |
912 | // estimated position in z of vertex | |
913 | zvtxTrks /= sumWeights; | |
914 | // estimated error | |
915 | ezvtxTrks = TMath::Sqrt(1./sumWeights); | |
916 | } else { | |
917 | zvtxTrks = zvtxSPD[ev]; | |
918 | } | |
919 | ||
920 | zvtx[ev] = zvtxTrks; | |
921 | ||
922 | } // loop on events | |
923 | ||
924 | ||
925 | // Write vertices to file | |
b2bca9d4 | 926 | WriteZvtx("zvtxTracks.dat",zvtx,n); |
056891c0 | 927 | delete [] zvtx; |
928 | delete [] zvtxSPD; | |
929 | ||
930 | gBenchmark->Stop(name); | |
931 | gBenchmark->Show(name); | |
932 | ||
933 | return 0; | |
934 | } | |
935 | //---------------------------------------------------------------------------- | |
936 | Bool_t VtxTrack(Double_t pt,Double_t d0rphi) { | |
937 | ||
938 | Double_t d0rphiRes = TMath::Sqrt(11.59*11.59+65.76*65.76/TMath::Power(pt,1.878))/10000.; | |
939 | Double_t nSigma = 3.; | |
940 | if(d0rphi<nSigma*d0rphiRes) { return kTRUE; } else { return kFALSE; } | |
941 | } | |
942 | //---------------------------------------------------------------------------- | |
943 | Double_t d0zRes(Double_t pt) { | |
944 | Double_t res = TMath::Sqrt(34.05*34.05+170.1*170.1/TMath::Power(pt,1.226)); | |
945 | return res/10000.; | |
946 | } | |
947 | //----------------------------------------------------------------------------- | |
948 | Int_t ITSFindTracksV2(const Char_t *galiceName,const Char_t *inName, | |
949 | const Char_t *inName2,const Char_t *outName, | |
b2bca9d4 | 950 | Bool_t *skipEvt,Option_t *vtxMode,Int_t n) { |
056891c0 | 951 | |
952 | ||
953 | cerr<<"\n*******************************************************************\n"; | |
954 | ||
955 | const Char_t *name="ITSFindTracksV2"; | |
956 | cerr<<'\n'<<name<<"...\n"; | |
957 | gBenchmark->Start(name); | |
958 | ||
959 | // Read vertices from file | |
960 | Double_t *zvtx = new Double_t[n]; | |
b2bca9d4 | 961 | Char_t *vtxfile="zvtxHeader.dat"; |
962 | ||
963 | const Char_t *header = strstr(vtxMode,"Header"); | |
964 | const Char_t *fastpp = strstr(vtxMode,"Fast"); | |
965 | const Char_t *spd = strstr(vtxMode,"SPD"); | |
966 | const Char_t *tracks = strstr(vtxMode,"Tracks"); | |
967 | ||
968 | if(header) vtxfile = "zvtxHeader.dat"; | |
969 | if(fastpp) vtxfile = "zvtxFastpp.dat"; | |
970 | if(spd) vtxfile = "zvtxSPD.dat"; | |
971 | if(tracks) vtxfile = "zvtxTracks.dat"; | |
056891c0 | 972 | |
b2bca9d4 | 973 | ReadZvtx(vtxfile,zvtx,n); |
056891c0 | 974 | |
975 | ||
976 | TFile *outFile = TFile::Open(outName,"recreate"); | |
977 | TFile *inFile = TFile::Open(inName); | |
978 | TFile *inFile2 = TFile::Open(inName2); | |
979 | ||
980 | AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom"); | |
981 | if(!geom) { cerr<<"can't get ITS geometry !\n"; return 1;} | |
982 | ||
983 | Int_t flag1stPass,flag2ndPass; | |
984 | ||
056891c0 | 985 | for(Int_t ev=0; ev<n; ev++){ |
056891c0 | 986 | if(skipEvt[ev]) continue; |
987 | cerr<<" --- Processing event "<<ev<<" ---"<<endl; | |
b2bca9d4 | 988 | AliITStrackerV2 tracker(geom,ev); |
056891c0 | 989 | |
990 | // set position of primary vertex | |
991 | Double_t vtx[3]; | |
992 | vtx[0]=0.; vtx[1]=0.; vtx[2]=zvtx[ev]; | |
993 | ||
994 | flag1stPass=1; // vtx constraint | |
995 | flag2ndPass=0; // no vtx constraint | |
996 | ||
997 | // no vtx constraint if vertex not found | |
998 | if(vtx[2]<-999.) { | |
999 | flag1stPass=0; | |
1000 | vtx[2]=0.; | |
1001 | } | |
1002 | ||
1003 | cerr<<"+++\n+++ Setting primary vertex z = "<<vtx[2]<< | |
1004 | " cm for ITS tracking\n+++\n"; | |
1005 | tracker.SetVertex(vtx); | |
88ce87da | 1006 | |
1007 | // setup vertex constraint in the two tracking passes | |
1008 | Int_t flags[2]; | |
056891c0 | 1009 | flags[0]=flag1stPass; |
88ce87da | 1010 | tracker.SetupFirstPass(flags); |
056891c0 | 1011 | flags[0]=flag2ndPass; |
88ce87da | 1012 | tracker.SetupSecondPass(flags); |
1013 | ||
056891c0 | 1014 | tracker.Clusters2Tracks(inFile,outFile); |
b2bca9d4 | 1015 | } |
056891c0 | 1016 | |
1017 | inFile->Close(); | |
1018 | inFile2->Close(); | |
1019 | outFile->Close(); | |
70521312 | 1020 | |
056891c0 | 1021 | delete [] zvtx; |
1022 | ||
88ce87da | 1023 | gBenchmark->Stop(name); |
1024 | gBenchmark->Show(name); | |
1025 | ||
056891c0 | 1026 | return 0; |
88ce87da | 1027 | } |
056891c0 | 1028 | //----------------------------------------------------------------------------- |
1029 | Int_t ITSMakeRefFile(const Char_t *galice,const Char_t *inname, | |
b2bca9d4 | 1030 | const Char_t *outname,Bool_t *skipEvt,Int_t n) { |
88ce87da | 1031 | |
1032 | ||
1033 | cerr<<"\n*******************************************************************\n"; | |
1034 | ||
1035 | Int_t rc=0; | |
1036 | const Char_t *name="ITSMakeRefFile"; | |
1037 | cerr<<'\n'<<name<<"...\n"; | |
1038 | gBenchmark->Start(name); | |
1039 | ||
1040 | ||
1041 | TFile *out = TFile::Open(outname,"recreate"); | |
1042 | TFile *trk = TFile::Open(inname); | |
1043 | TFile *kin = TFile::Open(galice); | |
1044 | ||
1045 | ||
1046 | // Get gAlice object from file | |
1047 | if(!(gAlice=(AliRun*)kin->Get("gAlice"))) { | |
1048 | cerr<<"gAlice has not been found on galice.root !\n"; | |
1049 | return 1; | |
1050 | } | |
1051 | ||
88ce87da | 1052 | Int_t label; |
1053 | TParticle *Part; | |
1054 | TParticle *Mum; | |
b2bca9d4 | 1055 | static RECTRACK rectrk; |
88ce87da | 1056 | |
1057 | ||
056891c0 | 1058 | for(Int_t ev=0; ev<n; ev++){ |
1059 | if(skipEvt[ev]) continue; | |
1060 | cerr<<" --- Processing event "<<ev<<" ---"<<endl; | |
88ce87da | 1061 | |
056891c0 | 1062 | gAlice->GetEvent(ev); |
88ce87da | 1063 | |
1064 | trk->cd(); | |
1065 | ||
1066 | // Tree with ITS tracks | |
1067 | char tname[100]; | |
056891c0 | 1068 | sprintf(tname,"TreeT_ITS_%d",ev); |
88ce87da | 1069 | |
1070 | TTree *tracktree=(TTree*)trk->Get(tname); | |
056891c0 | 1071 | if(!tracktree) continue; |
1072 | AliITStrackV2 *itstrack=new AliITStrackV2; | |
1073 | tracktree->SetBranchAddress("tracks",&itstrack); | |
88ce87da | 1074 | Int_t nentr=(Int_t)tracktree->GetEntries(); |
1075 | ||
1076 | // Tree for true track parameters | |
1077 | char ttname[100]; | |
056891c0 | 1078 | sprintf(ttname,"Tree_Ref_%d",ev); |
88ce87da | 1079 | TTree *reftree = new TTree(ttname,"Tree with true track params"); |
70521312 | 1080 | reftree->Branch("rectracks",&rectrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz"); |
88ce87da | 1081 | |
056891c0 | 1082 | for(Int_t i=0; i<nentr; i++) { |
88ce87da | 1083 | tracktree->GetEvent(i); |
1084 | label = TMath::Abs(itstrack->GetLabel()); | |
1085 | ||
1086 | Part = (TParticle*)gAlice->Particle(label); | |
1087 | rectrk.lab=label; | |
1088 | rectrk.pdg=Part->GetPdgCode(); | |
70521312 | 1089 | rectrk.mumlab = Part->GetFirstMother(); |
88ce87da | 1090 | if(Part->GetFirstMother()>=0) { |
1091 | Mum = (TParticle*)gAlice->Particle(Part->GetFirstMother()); | |
1092 | rectrk.mumpdg=Mum->GetPdgCode(); | |
1093 | } else { | |
1094 | rectrk.mumpdg=-1; | |
1095 | } | |
1096 | rectrk.Vx=Part->Vx(); | |
1097 | rectrk.Vy=Part->Vy(); | |
1098 | rectrk.Vz=Part->Vz(); | |
1099 | rectrk.Px=Part->Px(); | |
1100 | rectrk.Py=Part->Py(); | |
1101 | rectrk.Pz=Part->Pz(); | |
1102 | ||
1103 | reftree->Fill(); | |
056891c0 | 1104 | } // loop on tracks |
88ce87da | 1105 | |
1106 | out->cd(); | |
1107 | reftree->Write(); | |
1108 | ||
056891c0 | 1109 | delete itstrack; |
1110 | delete reftree; | |
1111 | } // loop on events | |
88ce87da | 1112 | |
1113 | trk->Close(); | |
1114 | kin->Close(); | |
1115 | out->Close(); | |
70521312 | 1116 | |
88ce87da | 1117 | gBenchmark->Stop(name); |
1118 | gBenchmark->Show(name); | |
1119 | ||
1120 | ||
1121 | return rc; | |
056891c0 | 1122 | } |
1123 | //----------------------------------------------------------------------------- | |
b2bca9d4 | 1124 | void WriteZvtx(const Char_t *name,Double_t *zvtx,Int_t n) { |
056891c0 | 1125 | |
b2bca9d4 | 1126 | FILE *f = fopen(name,"w"); |
1127 | for(Int_t i=0;i<n;i++) { | |
1128 | fprintf(f,"%f\n",zvtx[i]); | |
056891c0 | 1129 | } |
b2bca9d4 | 1130 | fclose(f); |
056891c0 | 1131 | |
1132 | return; | |
88ce87da | 1133 | } |
056891c0 | 1134 | //----------------------------------------------------------------------------- |
b2bca9d4 | 1135 | void ReadZvtx(const Char_t *name,Double_t *zvtx,Int_t n) { |
056891c0 | 1136 | |
b2bca9d4 | 1137 | FILE *f = fopen(name,"r"); |
1138 | for(Int_t i=0;i<n;i++) { | |
1139 | fscanf(f,"%lf",&zvtx[i]); | |
056891c0 | 1140 | } |
b2bca9d4 | 1141 | fclose(f); |
056891c0 | 1142 | |
1143 | return; | |
1144 | } | |
b2bca9d4 | 1145 | |
88ce87da | 1146 | |
1147 | ||
1148 | ||
1149 | ||
1150 | ||
1151 | ||
1152 | ||
1153 | ||
1154 |