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