1 void AliVertexerTracksTest(Double_t nSigma=3.,
2 Bool_t useMeanVtx=kFALSE,
8 TString outname="AliVertexerTracksTest.root",
9 Bool_t onlyfit=kFALSE) {
10 //------------------------------------------------------------------------
11 // Run vertex reconstruction and store results in histos and ntuple
12 //------------------------------------------------------------------------
13 // WITHOUT Kinematics.root
16 delete gAlice->GetRunLoader();
21 AliRunLoader *rl = AliRunLoader::Open("galice.root");
23 cerr<<"Can not open session"<<endl;
26 Int_t retval = rl->LoadgAlice();
28 cerr<<"LoadgAlice returned error"<<endl;
32 retval = rl->LoadHeader();
34 cerr<<"LoadHeader returned error"<<endl;
38 gAlice=rl->GetAliRun();
40 //rl->LoadKinematics();
42 // Get field from galice.root
43 AliMagF *fiel = (AliMagF*)gAlice->Field();
44 // Set the conversion constant between curvature and Pt
45 AliTracker::SetFieldMap(fiel,kTRUE);
47 //**** Switch on the PID class
50 Int_t nEvents = (Int_t)rl->GetNumberOfEvents();
55 Double_t diffX,diffY,diffZ;
56 Double_t diffXerr,diffYerr,diffZerr;
64 TH1F *hdiffX = new TH1F("hdiffX","x_{found} - x_{true} [#mum]",1000,-5000,5000);
65 TH1F *hdiffY = new TH1F("hdiffY","y_{found} - y_{true} [#mum]",1000,-5000,5000);
66 TH1F *hdiffZ = new TH1F("hdiffZ","z_{found} - z_{true} [#mum]",1000,-5000,5000);
68 TH1F *hdiffXerr = new TH1F("hdiffXerr","#Delta x/#sigma_{x}",100,-5,5);
69 TH1F *hdiffYerr = new TH1F("hdiffYerr","#Delta y/#sigma_{y}",100,-5,5);
70 TH1F *hdiffZerr = new TH1F("hdiffZerr","#Delta z/#sigma_{z}",100,-5,5);
73 TNtuple *ntVtxRes = new TNtuple("ntVtxRes","Vertex Resolution's Ntupla with multiplicity","triggered:ntracks:nitstracks:nitstracks5or6:diffX:diffY:diffZ:diffXerr:diffYerr:diffZerr:multiplicity:nc:zMC:chi2red");
77 // ----------- Create vertexer --------------------------
78 AliESDVertex *initVertex = 0;
80 // open the mean vertex
81 TFile *invtx = new TFile("AliESDVertexMean.root");
82 initVertex = (AliESDVertex*)invtx->Get("vtxmean");
86 Double_t pos[3]={0.5,0.5,0.};
87 Double_t err[3]={3.,3.,5.};
88 initVertex = new AliESDVertex(pos,err);
90 AliVertexerTracks *vertexer = new AliVertexerTracks();
91 vertexer->SetDebug(1); // set to 1 to see what it does
92 vertexer->SetVtxStart(initVertex);
93 if(!useMeanVtx) vertexer->SetNoConstraint();
94 if(onlyfit) vertexer->SetOnlyFitter();
95 vertexer->SetNSigmad0(nSigma);
96 if(!itsrefit) vertexer->SetITSrefitNotRequired();
97 if(!itsin) vertexer->SetITSNotRequired();
98 vertexer->SetMinTracks(mintracks);
99 vertexer->SetMinITSClusters(minitscls);
100 vertexer->SetMaxd0z0(maxd0z0);
101 //cout<<" Nsigma: "<<vertexer->GetNSigmad0()<<endl;
102 //----------------------------------------------------------
105 if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
106 printf("AliESDs.root not found!\n");
109 TFile *fin = TFile::Open("AliESDs.root");
110 TTree *esdTree = (TTree*)fin->Get("esdTree");
113 esdTree->SetBranchAddress("ESD",&esd);
114 Int_t events = esdTree->GetEntries();
115 printf("Number of events in ESD tree: %d\n",events);
119 for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
120 printf("-------------- EVENT %d --------------------\n",iev);
130 //=================================================
131 // LOOK IN THE SIMULATION
132 // true vertex position
133 Int_t npart = (Int_t)gAlice->GetEvent(iev);
134 printf("particles %d\n",npart);
135 AliHeader *header = (AliHeader*)rl->GetHeader();
136 AliGenPythiaEventHeader *pyheader = (AliGenPythiaEventHeader*)header->GenEventHeader();
137 pyheader->PrimaryVertex(o);
141 printf("True Vertex: (%f, %f, %f)\n",o[0],o[1],o[2]);
143 Double_t dNchdy = 0.;
145 Int_t nitstracks = 0;
146 Int_t nitstracks1 = 0;
147 Int_t nitstracks2 = 0;
148 Int_t nitstracks3 = 0;
149 Int_t nitstracks4 = 0;
150 Int_t nitstracks5 = 0;
151 Int_t nitstracks6 = 0;
152 Int_t nitstracks5or6 = 0;
155 esdTree->GetEvent(iev);
158 ULong64_t triggerMask = esd->GetTriggerMask();
159 if(triggerMask&32 && ((triggerMask&1) || (triggerMask&2))) triggered=1;
161 // get # of tracklets
162 AliMultiplicity *alimult = esd->GetMultiplicity();
164 ntracklets = alimult->GetNumberOfTracklets();
165 for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++)
166 if(alimult->GetDeltaPhi(l)<-9998.) ntracklets--;
170 printf("Number of SPD tracklets in ESD %d : %d\n",iev,ntracklets);
171 multiplicity = (Float_t)ntracklets;
173 Int_t ntracks = esd->GetNumberOfTracks();
174 printf("Number of tracks in ESD %d : %d\n",iev,ntracks);
176 ntVtxRes->Fill(triggered,ntracks,nitstracks,nitstracks5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,nc,truevtx[2],chi2red);
181 AliESDVertex *vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
183 nc = (Int_t)vertex->GetNContributors();
184 if(nc>=2) chi2red = vertex->GetChi2toNDF();
185 printf("nc = %d\n",nc);
189 for(Int_t itrk=0; itrk<ntracks; itrk++) {
190 AliESDtrack *esdtrack = (AliESDtrack*)esd->GetTrack(itrk);
191 UInt_t status = esdtrack->GetStatus();
192 // only tracks found also in ITS
193 if(! (status&AliESDtrack::kITSin) ) continue;
195 Int_t npts = (Int_t)esdtrack->GetNcls(0);
196 if(npts==6) {nitstracks6++;nitstracks5or6++;}
197 if(npts==5) {nitstracks5++;nitstracks5or6++;}
198 if(npts==4) nitstracks4++;
199 if(npts==3) nitstracks3++;
200 if(npts==2) nitstracks2++;
201 if(npts==1) nitstracks1++;
203 printf("Number of kITSin tracks in ESD %d : %d\n",iev,nitstracks);
204 printf(" 6 pts: %d\n",nitstracks6);
205 printf(" 5 pts: %d\n",nitstracks5);
206 printf(" 4 pts: %d\n",nitstracks4);
207 printf(" 3 pts: %d\n",nitstracks3);
208 printf(" 2 pts: %d\n",nitstracks2);
209 printf(" 1 pts: %d\n",nitstracks1);
214 vertex->GetSigmaXYZ(errvtx);
215 diffX = 10000.* (vtx[0] - truevtx[0]);
216 diffY = 10000.* (vtx[1] - truevtx[1]);
217 diffZ = 10000.* (vtx[2] - truevtx[2]);
221 diffXerr = diffX/errvtx[0]/10000.;
222 diffYerr = diffY/errvtx[1]/10000.;
223 diffZerr = diffZ/errvtx[2]/10000.;
224 hdiffXerr->Fill(diffXerr);
225 hdiffYerr->Fill(diffYerr);
226 hdiffZerr->Fill(diffZerr);
229 ntVtxRes->Fill(triggered,ntracks,nitstracks,nitstracks5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,(Float_t)nc,truevtx[2],chi2red);
231 } // END LOOP ON EVENTS
240 TFile *fout = new TFile(outname.Data(),"recreate");
252 //----------------------------------------------------------------------------
253 //------------------------------------------------------------------------
254 void VertexForOneEvent(Int_t iev=0,
256 Bool_t useMeanVtx=kFALSE) {
257 //------------------------------------------------------------------------
258 // Run vertex reconstruction for event iev
259 //------------------------------------------------------------------------
262 delete gAlice->GetRunLoader();
267 AliRunLoader *rl = AliRunLoader::Open("galice.root");
269 cerr<<"Can not open session"<<endl;
272 Int_t retval = rl->LoadgAlice();
274 cerr<<"LoadgAlice returned error"<<endl;
278 retval = rl->LoadHeader();
280 cerr<<"LoadHeader returned error"<<endl;
284 gAlice=rl->GetAliRun();
285 //rl->LoadKinematics();
286 // Get field from galice.root
287 AliMagF *fiel = (AliMagF*)gAlice->Field();
288 // Set the conversion constant between curvature and Pt
289 AliTracker::SetFieldMap(fiel,kTRUE);
294 Double_t diffX,diffY,diffZ;
295 Double_t diffXerr,diffYerr,diffZerr;
296 Double_t multiplicity;
299 // ----------- Create vertexer --------------------------
300 AliESDVertex *initVertex = 0;
302 // open the mean vertex
303 TFile *invtx = new TFile("AliESDVertexMean.root");
304 initVertex = (AliESDVertex*)invtx->Get("vtxmean");
308 Double_t pos[3]={0.,-0.,0.};
309 Double_t err[3]={3.,3.,5.};
310 initVertex = new AliESDVertex(pos,err);
312 AliVertexerTracks *vertexer = new AliVertexerTracks();
313 //vertexer->SetITSrefitNotRequired();
314 vertexer->SetVtxStart(initVertex);
315 if(!useMeanVtx) vertexer->SetNoConstraint();
316 //vertexer->SetMinTracks(2);
317 //vertexer->SetNSigmad0(nSigma);
318 //cout<<" Nsigma: "<<vertexer->GetNSigmad0()<<endl;
319 //vertexer->SetFinderAlgorithm(2);
320 //vertexer->SetDCAcut(0.1); // 1 mm
321 vertexer->SetDebug(1); // set to 1 to see what it does
322 //----------------------------------------------------------
324 if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
325 printf("AliESDs.root not found!\n");
328 TFile *fin = TFile::Open("AliESDs.root");
329 TTree *esdTree = (TTree*)fin->Get("esdTree");
332 esdTree->SetBranchAddress("ESD",&esd);
338 //=================================================
339 // LOOK IN THE SIMULATION
340 // true vertex position
341 Int_t npart = (Int_t)gAlice->GetEvent(iev);
342 printf("particles %d\n",npart);
343 AliHeader *header = (AliHeader*)rl->GetHeader();
344 AliGenPythiaEventHeader *pyheader = (AliGenPythiaEventHeader*)header->GenEventHeader();
345 pyheader->PrimaryVertex(o);
349 printf("True Vertex: (%f, %f, %f)\n",o[0],o[1],o[2]);
351 esdTree->GetEvent(iev);
352 Int_t ntracks = esd->GetNumberOfTracks();
353 printf("Number of tracks in ESD %d : %d\n",iev,ntracks);
356 AliESDVertex *vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
357 nc = (Int_t)vertex->GetNContributors();
360 Int_t nitstracks = 0;
361 Int_t nitstracks1 = 0;
362 Int_t nitstracks2 = 0;
363 Int_t nitstracks3 = 0;
364 Int_t nitstracks4 = 0;
365 Int_t nitstracks5 = 0;
366 Int_t nitstracks6 = 0;
367 Int_t nitstracks5or6 = 0;
370 for(Int_t itrk=0; itrk<ntracks; itrk++) {
371 AliESDtrack *esdtrack = (AliESDtrack*)esd->GetTrack(itrk);
372 UInt_t status = esdtrack->GetStatus();
373 // only tracks found also in ITS
374 if(! (status&AliESDtrack::kITSrefit) ) continue;
376 Int_t npts = (Int_t)esdtrack->GetNcls(0);
377 if(npts==6) {nitstracks6++;nitstracks5or6++;}
378 if(npts==5) {nitstracks5++;nitstracks5or6++;}
379 if(npts==4) nitstracks4++;
380 if(npts==3) nitstracks3++;
381 if(npts==2) nitstracks2++;
382 if(npts==1) nitstracks1++;
384 printf("Number of kITSrefit tracks in ESD %d : %d\n",iev,nitstracks);
385 printf(" 6 pts: %d\n",nitstracks6);
386 printf(" 5 pts: %d\n",nitstracks5);
387 printf(" 4 pts: %d\n",nitstracks4);
388 printf(" 3 pts: %d\n",nitstracks3);
389 printf(" 2 pts: %d\n",nitstracks2);
390 printf(" 1 pts: %d\n",nitstracks1);
400 //----------------------------------------------------------------------------
401 Double_t FitFunc(Double_t *x,Double_t *par) {
402 return par[0]+par[1]/TMath::Sqrt(x[0]);
404 Int_t GetBin(Float_t mult) {
405 if(mult<5.) return 0;
406 if(mult<7.) return 1;
407 if(mult<12.) return 2;
408 if(mult<15.) return 3;
409 if(mult<22.) return 4;
412 Int_t GetBinTracklets(Float_t tracklets) {
413 if(tracklets<2.*5.) return 0;
414 if(tracklets<2.*7.) return 1;
415 if(tracklets<2.*12.) return 2;
416 if(tracklets<2.*15.) return 3;
417 if(tracklets<2.*22.) return 4;
420 Int_t GetBinZ(Float_t z) {
437 //----------------------------------------------------------------------------
438 void PlotVtxRes(TString inname="AliVertexerTracksTest.root") {
439 //---------------------------------------------------------------------
440 // Plot efficiency, resolutions, pulls
441 // [at least 10k pp events are needed]
442 //---------------------------------------------------------------------
446 Float_t ncminforzshift=1.;
448 gStyle->SetOptStat(0);
449 gStyle->SetOptFit(10001);
454 Float_t nESDTrks,nTrks,nTrks5or6,ntracklets;
455 Float_t diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,zMC,nc;
462 TH1F *hdx = new TH1F("hdx","",50,-600,600);
463 hdx->SetXTitle("#Delta x [#mu m]");
464 hdx->SetYTitle("events");
466 TH1F *hdy = new TH1F("hdy","",50,-600,600);
467 hdy->SetXTitle("#Delta y [#mu m]");
468 hdy->SetYTitle("events");
470 TH1F *hdz = new TH1F("hdz","",200,-600,600);
471 hdz->SetXTitle("#Delta z [#mu m]");
472 hdz->SetYTitle("events");
474 TH1F *hzMCVtx = new TH1F("hzMCVtx","events w/ vertex",30,-15,15);
475 hzMCVtx->SetXTitle("z_{true} [cm]");
476 hzMCVtx->SetYTitle("events");
479 TH1F *hzMCNoVtx = new TH1F("hzMCNoVtx","events w/o vertex",30,-15,15);
480 hzMCNoVtx->SetXTitle("z_{true} [cm]");
481 hzMCNoVtx->SetYTitle("events");
484 TH1F *hTrackletsVtx = new TH1F("hTrackletsVtx","events w/ vertex",51,-0.5,50.5);
485 hTrackletsVtx->SetXTitle("SPD tracklets");
486 hTrackletsVtx->SetYTitle("events");
488 TH1F *hTrackletsNoVtx = new TH1F("hTrackletsNoVtx","events w/o vertex",51,-0.5,50.5);
489 hTrackletsNoVtx->SetXTitle("SPD tracklets");
490 hTrackletsNoVtx->SetYTitle("events");
492 TH1F *hTracksVtx = new TH1F("hTracksVtx","events w/ vertex",51,-0.5,50.5);
493 hTracksVtx->SetXTitle("Number of reco tracks in TPC&ITS");
494 hTracksVtx->SetYTitle("events");
496 TH1F *hTracksNoVtx = new TH1F("hTracksNoVtx","events w/o vertex",51,-0.5,50.5);
497 hTracksNoVtx->SetXTitle("Number of reco tracks in TPC&ITS");
498 hTracksNoVtx->SetYTitle("events");
500 TH1F *hTracksVtx5or6 = new TH1F("hTracksVtx5or6","events w/ vertex",51,-0.5,50.5);
501 hTracksVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(cls>=5)");
502 hTracksVtx5or6->SetYTitle("events");
504 TH1F *hTracksNoVtx5or6 = new TH1F("hTracksNoVtx5or6","events w/o vertex",51,-0.5,50.5);
505 hTracksNoVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(cls>=5)");
506 hTracksNoVtx5or6->SetYTitle("events");
508 TH1F *hPullsx = new TH1F("hPullsx","",50,-7.,7.);
509 hPullsx->SetXTitle("#Delta x/#sqrt{<x,x>}");
510 hPullsx->SetYTitle("events");
512 TH1F *hPullsy = new TH1F("hPullsy","",50,-7.,7.);
513 hPullsy->SetXTitle("#Delta y/#sqrt{<y,y>}");
514 hPullsy->SetYTitle("events");
516 TH1F *hPullsz = new TH1F("hPullsz","",50,-7.,7.);
517 hPullsz->SetXTitle("#Delta z/#sqrt{<z,z>}");
518 hPullsz->SetYTitle("events");
520 TProfile *hntrackletsVSzMC = new TProfile("hntrackletsVSzMC","",30,-15,15,0,30);
521 hntrackletsVSzMC->SetXTitle("z_{true} [cm]");
522 hntrackletsVSzMC->SetYTitle("<SPD tracklets>");
524 TProfile *hnitstracksVSzMC = new TProfile("hnitstracksVSzMC","",30,-15,15,0,30);
525 hnitstracksVSzMC->SetXTitle("z_{true} [cm]");
526 hnitstracksVSzMC->SetYTitle("<tracks TPC&ITS>");
528 TProfile *hnitstracks5or6VSzMC = new TProfile("hnitstracks5or6VSzMC","",30,-15,15,0,30);
529 hnitstracks5or6VSzMC->SetXTitle("z_{true} [cm]");
530 hnitstracks5or6VSzMC->SetYTitle("<tracks TPC&ITS(cls>=5)>");
532 Double_t mult[6]={0,0,0,0,0,0};
533 Double_t mult2[6]={0,0,0,0,0,0};
534 Double_t sigmamult[6]={0,0,0,0,0,0};
535 Double_t sigmamult2[6]={0,0,0,0,0,0};
536 Double_t tottracks[6]={0,0,0,0,0,0};
537 Double_t fittrks[6]={0,0,0,0,0,0};
538 Double_t tracks[6]={0,0,0,0,0,0};
539 Double_t etracks[6]={0,0,0,0,0,0};
540 Double_t xres[6]={0,0,0,0,0,0};
541 Double_t exres[6]={0,0,0,0,0,0};
542 Double_t yres[6]={0,0,0,0,0,0};
543 Double_t eyres[6]={0,0,0,0,0,0};
544 Double_t zres[6]={0,0,0,0,0,0};
545 Double_t ezres[6]={0,0,0,0,0,0};
546 Double_t xpull[6]={0,0,0,0,0,0};
547 Double_t expull[6]={0,0,0,0,0,0};
548 Double_t ypull[6]={0,0,0,0,0,0};
549 Double_t eypull[6]={0,0,0,0,0,0};
550 Double_t zpull[6]={0,0,0,0,0,0};
551 Double_t ezpull[6]={0,0,0,0,0,0};
552 Int_t evts[6]={0,0,0,0,0,0};
553 Int_t totevts[6]={0,0,0,0,0,0};
554 Int_t vtx[6]={0,0,0,0,0,0};
555 Int_t all[6]={0,0,0,0,0,0};
556 Double_t probVtx[6]={0,0,0,0,0,0};
557 Double_t eprobVtx[6]={0,0,0,0,0,0};
562 TH1F *hdx0 = new TH1F("hdx0","x resolution - bin 0",50,-500,500);
563 TH1F *hdx1 = new TH1F("hdx1","x resolution - bin 1",50,-500,500);
564 TH1F *hdx2 = new TH1F("hdx2","x resolution - bin 2",50,-500,500);
565 TH1F *hdx3 = new TH1F("hdx3","x resolution - bin 3",50,-400,400);
566 TH1F *hdx4 = new TH1F("hdx4","x resolution - bin 4",50,-300,300);
567 TH1F *hdx5 = new TH1F("hdx5","x resolution - bin 5",50,-300,300);
568 TH1F *hdy0 = new TH1F("hdy0","y resolution - bin 0",50,-500,500);
569 TH1F *hdy1 = new TH1F("hdy1","y resolution - bin 1",50,-500,500);
570 TH1F *hdy2 = new TH1F("hdy2","y resolution - bin 2",50,-500,500);
571 TH1F *hdy3 = new TH1F("hdy3","y resolution - bin 3",50,-400,400);
572 TH1F *hdy4 = new TH1F("hdy4","y resolution - bin 4",50,-300,300);
573 TH1F *hdy5 = new TH1F("hdy5","y resolution - bin 5",50,-300,300);
574 TH1F *hdz0 = new TH1F("hdz0","z resolution - bin 0",50,-1000,1000);
575 TH1F *hdz1 = new TH1F("hdz1","z resolution - bin 1",50,-800,800);
576 TH1F *hdz2 = new TH1F("hdz2","z resolution - bin 2",50,-800,800);
577 TH1F *hdz3 = new TH1F("hdz3","z resolution - bin 3",50,-600,600);
578 TH1F *hdz4 = new TH1F("hdz4","z resolution - bin 4",50,-500,500);
579 TH1F *hdz5 = new TH1F("hdz5","z resolution - bin 5",50,-500,500);
581 TH1F *hdz_z = NULL; hdz_z = new TH1F[15];
582 for(Int_t l=0;l<15;l++) hdz_z[l] = *hdz;
585 TH1F *hpx0 = new TH1F("hpx0","x pull - bin 0",50,-7,7);
586 TH1F *hpx1 = new TH1F("hpx1","x pull - bin 1",50,-7,7);
587 TH1F *hpx2 = new TH1F("hpx2","x pull - bin 2",50,-7,7);
588 TH1F *hpx3 = new TH1F("hpx3","x pull - bin 3",50,-7,7);
589 TH1F *hpx4 = new TH1F("hpx4","x pull - bin 4",50,-7,7);
590 TH1F *hpx5 = new TH1F("hpx5","x pull - bin 5",50,-7,7);
591 TH1F *hpy0 = new TH1F("hpy0","y pull - bin 0",50,-7,7);
592 TH1F *hpy1 = new TH1F("hpy1","y pull - bin 1",50,-7,7);
593 TH1F *hpy2 = new TH1F("hpy2","y pull - bin 2",50,-7,7);
594 TH1F *hpy3 = new TH1F("hpy3","y pull - bin 3",50,-7,7);
595 TH1F *hpy4 = new TH1F("hpy4","y pull - bin 4",50,-7,7);
596 TH1F *hpy5 = new TH1F("hpy5","y pull - bin 5",50,-7,7);
597 TH1F *hpz0 = new TH1F("hpz0","z pull - bin 0",50,-7,7);
598 TH1F *hpz1 = new TH1F("hpz1","z pull - bin 1",50,-7,7);
599 TH1F *hpz2 = new TH1F("hpz2","z pull - bin 2",50,-7,7);
600 TH1F *hpz3 = new TH1F("hpz3","z pull - bin 3",50,-7,7);
601 TH1F *hpz4 = new TH1F("hpz4","z pull - bin 4",50,-7,7);
602 TH1F *hpz5 = new TH1F("hpz5","z pull - bin 5",50,-7,7);
605 TH1F *hmult = new TH1F("hmult","hmult",50,-0.5,49.5);
607 TFile *in = new TFile(inname.Data());
608 TNtuple *nt = (TNtuple*)in->Get("ntVtxRes");
609 nt->SetBranchAddress("triggered",&triggered);
610 nt->SetBranchAddress("ntracks",&nESDTrks);
611 nt->SetBranchAddress("nitstracks",&nTrks);
612 nt->SetBranchAddress("nitstracks5or6",&nTrks5or6);
613 nt->SetBranchAddress("diffX",&diffX);
614 nt->SetBranchAddress("diffY",&diffY);
615 nt->SetBranchAddress("diffZ",&diffZ);
616 nt->SetBranchAddress("diffXerr",&diffXerr);
617 nt->SetBranchAddress("diffYerr",&diffYerr);
618 nt->SetBranchAddress("diffZerr",&diffZerr);
619 nt->SetBranchAddress("multiplicity",&ntracklets);
620 nt->SetBranchAddress("zMC",&zMC);
621 nt->SetBranchAddress("nc",&nc);
622 Int_t entries = (Int_t)nt->GetEntries();
625 for(Int_t i=0; i<entries; i++) {
626 nbytes += nt->GetEvent(i);
627 // only triggered events
628 if(triggered<0.5) continue;
629 // consider only events with zMCmin<|zMC|<zMCmax
630 if(TMath::Abs(zMC)<zMCmin || TMath::Abs(zMC)>zMCmax) continue;
632 bin = GetBinTracklets(ntracklets);
634 mult[bin] += ntracklets;
635 hmult->Fill(ntracklets);
637 hntrackletsVSzMC->Fill(zMC,ntracklets);
638 hnitstracksVSzMC->Fill(zMC,nTrks);
639 hnitstracks5or6VSzMC->Fill(zMC,nTrks5or6);
640 if(TMath::Abs(diffX) < 1000.) { // vtx found - closer than 1 mm
642 mult2[bin] += ntracklets;
644 tottracks[bin] += nTrks;
646 if(bin==0) hdx0->Fill(diffX);
647 if(bin==1) hdx1->Fill(diffX);
648 if(bin==2) hdx2->Fill(diffX);
649 if(bin==3) hdx3->Fill(diffX);
650 if(bin==4) hdx4->Fill(diffX);
651 if(bin==5) hdx5->Fill(diffX);
652 if(bin==0) hdy0->Fill(diffY);
653 if(bin==1) hdy1->Fill(diffY);
654 if(bin==2) hdy2->Fill(diffY);
655 if(bin==3) hdy3->Fill(diffY);
656 if(bin==4) hdy4->Fill(diffY);
657 if(bin==5) hdy5->Fill(diffY);
658 if(bin==0) hdz0->Fill(diffZ);
659 if(bin==1) hdz1->Fill(diffZ);
660 if(bin==2) hdz2->Fill(diffZ);
661 if(bin==3) hdz3->Fill(diffZ);
662 if(bin==4) hdz4->Fill(diffZ);
663 if(bin==5) hdz5->Fill(diffZ);
667 if(nc>=ncminforzshift) hdz_z[zbin].Fill(diffZ);
669 hTrackletsVtx->Fill(ntracklets);
670 hTracksVtx->Fill(nTrks);
671 hTracksVtx5or6->Fill(nTrks5or6);
672 hPullsx->Fill(diffXerr);
673 hPullsy->Fill(diffYerr);
674 hPullsz->Fill(diffZerr);
675 if(bin==0) hpx0->Fill(diffXerr);
676 if(bin==1) hpx1->Fill(diffXerr);
677 if(bin==2) hpx2->Fill(diffXerr);
678 if(bin==3) hpx3->Fill(diffXerr);
679 if(bin==4) hpx4->Fill(diffXerr);
680 if(bin==5) hpx5->Fill(diffXerr);
681 if(bin==0) hpy0->Fill(diffYerr);
682 if(bin==1) hpy1->Fill(diffYerr);
683 if(bin==2) hpy2->Fill(diffYerr);
684 if(bin==3) hpy3->Fill(diffYerr);
685 if(bin==4) hpy4->Fill(diffYerr);
686 if(bin==5) hpy5->Fill(diffYerr);
687 if(bin==0) hpz0->Fill(diffZerr);
688 if(bin==1) hpz1->Fill(diffZerr);
689 if(bin==2) hpz2->Fill(diffZerr);
690 if(bin==3) hpz3->Fill(diffZerr);
691 if(bin==4) hpz4->Fill(diffZerr);
692 if(bin==5) hpz5->Fill(diffZerr);
695 hzMCNoVtx->Fill(zMC);
696 hTrackletsNoVtx->Fill(ntracklets);
697 hTracksNoVtx->Fill(nTrks);
698 hTracksNoVtx5or6->Fill(nTrks5or6);
704 for(Int_t k=0;k<6;k++) {
706 mult[k] /= totevts[k];
707 tottracks[k] /= evts[k];
708 probVtx[k] = (Double_t)vtx[k]/all[k];
709 eprobVtx[k] = (Double_t)vtx[k]/all[k]/all[k]+(Double_t)vtx[k]*vtx[k]/all[k]/all[k]/all[k];
710 eprobVtx[k] = TMath::Sqrt(eprobVtx[k]);
713 // calculate spread in multiplicity
714 for(Int_t i=0; i<entries; i++) {
715 nbytes += nt->GetEvent(i);
716 bin = GetBinTracklets(ntracklets);
717 sigmamult[bin] += (ntracklets-mult[bin])*(ntracklets-mult[bin])/totevts[bin];
718 if(diffX < 90000.) sigmamult2[bin] += (ntracklets-mult2[bin])*(ntracklets-mult2[bin])/evts[bin];
721 for(Int_t k=0;k<6;k++) {
722 sigmamult[k] = TMath::Sqrt(sigmamult[k]);
723 sigmamult2[k] = TMath::Sqrt(sigmamult2[k]);
726 // fraction of found vertices
727 Float_t frac = (Float_t)nEvVtx/(nEvVtx+nEvNoVtx);
728 printf(" Fraction of events with vertex: %f\n",frac);
732 TF1 *g = new TF1("g","gaus");
734 TF1 *line = new TF1("line","pol1");
735 line->SetLineWidth(3);
736 TF1 *func = new TF1("func",FitFunc,2,80,2);
737 func->SetLineWidth(1);
739 Double_t x1,y1,sigmafit;
745 gStyle->SetOptFit(1111);
748 TCanvas *c0 = new TCanvas("c0","c0",0,0,1000,400);
751 hntrackletsVSzMC->SetMinimum(0);
752 hntrackletsVSzMC->SetMaximum(14);
753 hntrackletsVSzMC->Draw("profs");
755 hnitstracksVSzMC->SetMinimum(0);
756 hnitstracksVSzMC->SetMaximum(14);
757 hnitstracksVSzMC->Draw("profs");
759 hnitstracks5or6VSzMC->SetMinimum(0);
760 hnitstracks5or6VSzMC->SetMaximum(14);
761 hnitstracks5or6VSzMC->Draw("profs");
764 // Tracks in ITS for events w/ and w/o vertex
765 TCanvas *c1 = new TCanvas("c1","c1",0,0,1000,500);
770 hTracksNoVtx->Draw();
772 // probability vs ntracklets
773 TCanvas *c1a = new TCanvas("c1a","c1a",0,0,500,500);
774 TH1F *hTotTracklets = (TH1F*)hTrackletsVtx->Clone("hTotTracklets");
775 hTotTracklets->Add(hTrackletsNoVtx);
776 TH1F *hProbTracklets = (TH1F*)hTrackletsVtx->Clone("hProbTracklets");
777 hProbTracklets->Divide(hTotTracklets);
778 hProbTracklets->SetYTitle("Probability");
779 hProbTracklets->SetTitle("Probability to find the vertex");
780 hProbTracklets->Draw();
782 // probability vs ntracks
783 TCanvas *c1b = new TCanvas("c1b","c1b",0,0,500,500);
784 TH1F *hTotTracks = (TH1F*)hTracksVtx->Clone("hTotTracks");
785 hTotTracks->Add(hTracksNoVtx);
786 TH1F *hProbTracks = (TH1F*)hTracksVtx->Clone("hProbTracks");
787 hProbTracks->Divide(hTotTracks);
788 hProbTracks->SetYTitle("Probability");
789 hProbTracks->SetTitle("Probability to find the vertex");
792 TCanvas *c1c = new TCanvas("c1c","c1c",0,0,500,500);
793 TH1F *hTotTracks5or6 = (TH1F*)hTracksVtx5or6->Clone("hTotTracks5or6");
794 hTotTracks5or6->Add(hTracksNoVtx5or6);
795 TH1F *hProbTracks5or6 = (TH1F*)hTracksVtx5or6->Clone("hProbTracks5or6");
796 hProbTracks5or6->Divide(hTotTracks5or6);
797 hProbTracks5or6->SetYTitle("Probability");
798 hProbTracks5or6->SetTitle("Probability to find the vertex");
799 hProbTracks5or6->Draw();
802 // probability vs zMC
803 TCanvas *c1e = new TCanvas("c1e","c1e",0,0,500,500);
804 TH1F *hTotzMC = (TH1F*)hzMCVtx->Clone("hTotzMC");
805 hTotzMC->Add(hzMCNoVtx);
806 TH1F *hProbzMC = (TH1F*)hzMCVtx->Clone("hProbzMC");
807 hProbzMC->Divide(hTotzMC);
808 hProbzMC->SetYTitle("Probability");
809 hProbzMC->SetTitle("Probability to find the vertex");
812 // Global resolutions
813 TCanvas *c2 = new TCanvas("c2","c2",0,0,1000,400);
817 g->SetRange(-1.*hdx->GetRMS(),+1.*hdx->GetRMS());
821 g->SetRange(-1.*hdy->GetRMS(),+1.*hdy->GetRMS());
825 g->SetRange(-1.*hdz->GetRMS(),+1.*hdz->GetRMS());
829 TCanvas *c4 = new TCanvas("c4","c4",0,0,1000,400);
833 g->SetRange(-3.*hPullsx->GetRMS(),+3.*hPullsx->GetRMS());
834 hPullsx->Fit("g","R");
837 g->SetRange(-3.*hPullsy->GetRMS(),+3.*hPullsy->GetRMS());
838 hPullsy->Fit("g","R");
841 g->SetRange(-3.*hPullsz->GetRMS(),+3.*hPullsz->GetRMS());
842 hPullsz->Fit("g","R");
844 // probability VS multiplicity
845 TCanvas *c5 = new TCanvas("c5","c5",0,0,600,500);
846 TH2F *fr5 = new TH2F("fr5","",2,0,80,2,0,1.1);
847 fr5->SetXTitle("SPD tracklets");
848 fr5->SetYTitle("Probability");
850 TGraphErrors *gr5 = new TGraphErrors(6,mult,probVtx,sigmamult,eprobVtx);
852 gr5->SetMarkerStyle(21);
855 // resolutions VS multiplicity
858 TCanvas *ax = new TCanvas("ax","ax",0,0,900,700);
862 g->SetRange(-3.*hdx0->GetRMS(),+3.*hdx0->GetRMS());
864 xres[0]=g->GetParameter(2);
865 exres[0]=g->GetParError(2);
868 g->SetRange(-3.*hdx1->GetRMS(),+3.*hdx1->GetRMS());
870 xres[1]=g->GetParameter(2);
871 exres[1]=g->GetParError(2);
874 g->SetRange(-3.*hdx2->GetRMS(),+3.*hdx2->GetRMS());
876 xres[2]=g->GetParameter(2);
877 exres[2]=g->GetParError(2);
880 g->SetRange(-3.*hdx3->GetRMS(),+3.*hdx3->GetRMS());
882 xres[3]=g->GetParameter(2);
883 exres[3]=g->GetParError(2);
886 g->SetRange(-3.*hdx4->GetRMS(),+3.*hdx4->GetRMS());
888 xres[4]=g->GetParameter(2);
889 exres[4]=g->GetParError(2);
892 g->SetRange(-3.*hdx5->GetRMS(),+3.*hdx5->GetRMS());
894 xres[5]=g->GetParameter(2);
895 exres[5]=g->GetParError(2);
897 TCanvas *ay = new TCanvas("ay","ay",0,0,900,700);
901 g->SetRange(-3.*hdy0->GetRMS(),+3.*hdy0->GetRMS());
903 yres[0]=g->GetParameter(2);
904 eyres[0]=g->GetParError(2);
907 g->SetRange(-3.*hdy1->GetRMS(),+3.*hdy1->GetRMS());
909 yres[1]=g->GetParameter(2);
910 eyres[1]=g->GetParError(2);
913 g->SetRange(-3.*hdy2->GetRMS(),+3.*hdy2->GetRMS());
915 yres[2]=g->GetParameter(2);
916 eyres[2]=g->GetParError(2);
919 g->SetRange(-3.*hdy3->GetRMS(),+3.*hdy3->GetRMS());
921 yres[3]=g->GetParameter(2);
922 eyres[3]=g->GetParError(2);
925 g->SetRange(-3.*hdy4->GetRMS(),+3.*hdy4->GetRMS());
927 yres[4]=g->GetParameter(2);
928 eyres[4]=g->GetParError(2);
931 g->SetRange(-3.*hdy5->GetRMS(),+3.*hdy5->GetRMS());
933 yres[5]=g->GetParameter(2);
934 eyres[5]=g->GetParError(2);
936 TCanvas *az = new TCanvas("az","az",0,0,900,700);
940 g->SetRange(-3.*hdz0->GetRMS(),+3.*hdz0->GetRMS());
942 zres[0]=g->GetParameter(2);
943 ezres[0]=g->GetParError(2);
946 g->SetRange(-3.*hdz1->GetRMS(),+3.*hdz1->GetRMS());
948 zres[1]=g->GetParameter(2);
949 ezres[1]=g->GetParError(2);
952 g->SetRange(-3.*hdz2->GetRMS(),+3.*hdz2->GetRMS());
954 zres[2]=g->GetParameter(2);
955 ezres[2]=g->GetParError(2);
958 g->SetRange(-3.*hdz3->GetRMS(),+3.*hdz3->GetRMS());
960 zres[3]=g->GetParameter(2);
961 ezres[3]=g->GetParError(2);
964 g->SetRange(-3.*hdz4->GetRMS(),+3.*hdz4->GetRMS());
966 zres[4]=g->GetParameter(2);
967 ezres[4]=g->GetParError(2);
970 g->SetRange(-3.*hdz5->GetRMS(),+3.*hdz5->GetRMS());
972 zres[5]=g->GetParameter(2);
973 ezres[5]=g->GetParError(2);
977 // pulls VS multiplicity
980 TCanvas *bx = new TCanvas("bx","bx",0,0,900,700);
984 g->SetRange(-3.*hpx0->GetRMS(),+3.*hpx0->GetRMS());
986 xpull[0]=g->GetParameter(2);
987 expull[0]=g->GetParError(2);
990 g->SetRange(-3.*hpx1->GetRMS(),+3.*hpx1->GetRMS());
992 xpull[1]=g->GetParameter(2);
993 expull[1]=g->GetParError(2);
996 g->SetRange(-3.*hpx2->GetRMS(),+3.*hpx2->GetRMS());
998 xpull[2]=g->GetParameter(2);
999 expull[2]=g->GetParError(2);
1002 g->SetRange(-3.*hpx3->GetRMS(),+3.*hpx3->GetRMS());
1004 xpull[3]=g->GetParameter(2);
1005 expull[3]=g->GetParError(2);
1008 g->SetRange(-3.*hpx4->GetRMS(),+3.*hpx4->GetRMS());
1010 xpull[4]=g->GetParameter(2);
1011 expull[4]=g->GetParError(2);
1014 g->SetRange(-3.*hpx5->GetRMS(),+3.*hpx5->GetRMS());
1016 xpull[5]=g->GetParameter(2);
1017 expull[5]=g->GetParError(2);
1019 TCanvas *by = new TCanvas("by","by",0,0,900,700);
1023 g->SetRange(-3.*hpy0->GetRMS(),+3.*hpy0->GetRMS());
1025 ypull[0]=g->GetParameter(2);
1026 eypull[0]=g->GetParError(2);
1029 g->SetRange(-3.*hpy1->GetRMS(),+3.*hpy1->GetRMS());
1031 ypull[1]=g->GetParameter(2);
1032 eypull[1]=g->GetParError(2);
1035 g->SetRange(-3.*hpy2->GetRMS(),+3.*hpy2->GetRMS());
1037 ypull[2]=g->GetParameter(2);
1038 eypull[2]=g->GetParError(2);
1041 g->SetRange(-3.*hpy3->GetRMS(),+3.*hpy3->GetRMS());
1043 ypull[3]=g->GetParameter(2);
1044 eypull[3]=g->GetParError(2);
1047 g->SetRange(-3.*hpy4->GetRMS(),+3.*hpy4->GetRMS());
1049 ypull[4]=g->GetParameter(2);
1050 eypull[4]=g->GetParError(2);
1053 g->SetRange(-3.*hpy5->GetRMS(),+3.*hpy5->GetRMS());
1055 ypull[5]=g->GetParameter(2);
1056 eypull[5]=g->GetParError(2);
1058 TCanvas *bz = new TCanvas("bz","bz",0,0,900,700);
1062 g->SetRange(-3.*hpz0->GetRMS(),+3.*hpz0->GetRMS());
1064 zpull[0]=g->GetParameter(2);
1065 ezpull[0]=g->GetParError(2);
1068 g->SetRange(-3.*hpz1->GetRMS(),+3.*hpz1->GetRMS());
1070 zpull[1]=g->GetParameter(2);
1071 ezpull[1]=g->GetParError(2);
1074 g->SetRange(-3.*hpz2->GetRMS(),+3.*hpz2->GetRMS());
1076 zpull[2]=g->GetParameter(2);
1077 ezpull[2]=g->GetParError(2);
1080 g->SetRange(-3.*hpz3->GetRMS(),+3.*hpz3->GetRMS());
1082 zpull[3]=g->GetParameter(2);
1083 ezpull[3]=g->GetParError(2);
1086 g->SetRange(-3.*hpz4->GetRMS(),+3.*hpz4->GetRMS());
1088 zpull[4]=g->GetParameter(2);
1089 ezpull[4]=g->GetParError(2);
1092 g->SetRange(-3.*hpz5->GetRMS(),+3.*hpz5->GetRMS());
1094 zpull[5]=g->GetParameter(2);
1095 ezpull[5]=g->GetParError(2);
1097 gStyle->SetOptFit(0);
1099 // res VS ntracklets
1100 TCanvas *c6 = new TCanvas("c6","c6",0,0,600,500);
1101 TH2F *fr6 = new TH2F("fr6","",2,0,80,2,0,240);
1102 fr6->SetXTitle("SPD tracklets");
1103 fr6->SetYTitle("#sigma [#mu m]");
1105 sigmamult2[0]=sigmamult2[1]=sigmamult2[2]=sigmamult2[3]=sigmamult2[4]=sigmamult2[5]=0.;
1106 TGraphErrors *gr6x = new TGraphErrors(6,mult2,xres,sigmamult2,exres);
1108 gr6x->SetMarkerStyle(22);
1109 //gr6x->Fit("func","E,R");
1110 TGraphErrors *gr6z = new TGraphErrors(6,mult2,zres,sigmamult2,ezres);
1112 gr6z->SetMarkerStyle(26);
1113 //gr6z->Fit("func","E,R");
1114 TLegend *leg6 = new TLegend(0.6,0.75,0.9,0.9);
1115 leg6->AddEntry(gr6x,"x/y coordinate","p");
1116 leg6->AddEntry(gr6z,"z coordinate","p");
1117 leg6->SetFillStyle(0);
1118 leg6->SetBorderSize(0);
1121 // pull VS ntracklets
1122 TCanvas *c8 = new TCanvas("c8","c8",0,0,600,500);
1123 TH2F *fr8 = new TH2F("fr8","",2,0,80,2,0,2.);
1124 fr8->SetXTitle("SPD tracklets");
1125 fr8->SetYTitle("pull");
1127 TGraphErrors *gr8x = new TGraphErrors(6,mult2,xpull,sigmamult2,expull);
1129 gr8x->SetMarkerStyle(22);
1130 TGraphErrors *gr8z = new TGraphErrors(6,mult2,zpull,sigmamult2,ezpull);
1132 gr8z->SetMarkerStyle(26);
1133 TLegend *leg8 = new TLegend(0.6,0.75,0.9,0.9);
1134 leg8->AddEntry(gr6x,"x/y coordinate","p");
1135 leg8->AddEntry(gr6z,"z coordinate","p");
1136 leg8->SetFillStyle(0);
1137 leg8->SetBorderSize(0);
1139 TLine *l8 = new TLine(0,1,80,1);
1140 l8->SetLineStyle(2);
1144 // mean z res VS zMC
1145 Float_t zmc[15]={-14,-12,-10,-8,-6,-4,-2,0,2,4,6,8,10,12,14};
1147 Float_t dzmean[15],edzmean[15];
1149 TCanvas *zz = new TCanvas("zz","zz",0,0,1000,1000);
1154 g->SetRange(-3.*hdz_z[l].GetRMS(),+3.*hdz_z[l].GetRMS());
1155 hdz_z[l].Fit("g","Q");
1156 dzmean[l]=g->GetParameter(1);
1157 edzmean[l]=g->GetParError(1);
1160 TCanvas *zzz = new TCanvas("zzz","zzz",0,0,600,500);
1161 TH2F *fr9 = new TH2F("fr9","",2,-15,15,2,-50,50);
1162 fr9->SetXTitle("z_{true} [cm]");
1163 fr9->SetYTitle("<z_{rec} - z_{true}> [#mu m]");
1165 TGraphErrors *grz = new TGraphErrors(15,zmc,dzmean,ezmc,edzmean);
1167 grz->SetMarkerStyle(20);
1172 //----------------------------------------------------------------------------
1173 void SaveFigures(TString name="dummy") {
1177 namefull = name.Data();
1178 namefull.Prepend("plots/residualsx_");
1179 namefull.Append(".gif");
1180 ax->Print(namefull.Data());
1181 namefull = name.Data();
1182 namefull.Prepend("plots/residualsx_");
1183 namefull.Append(".ps");
1184 ax->Print(namefull.Data());
1185 namefull = name.Data();
1186 namefull.Prepend("plots/residualsx_");
1187 namefull.Append(".eps");
1188 ax->Print(namefull.Data());
1189 namefull = name.Data();
1190 namefull.Prepend("plots/residualsx_");
1191 namefull.Append(".root");
1192 ax->Print(namefull.Data());
1194 namefull = name.Data();
1195 namefull.Prepend("plots/pullsx_");
1196 namefull.Append(".gif");
1197 bx->Print(namefull.Data());
1198 namefull = name.Data();
1199 namefull.Prepend("plots/pullsx_");
1200 namefull.Append(".ps");
1201 bx->Print(namefull.Data());
1202 namefull = name.Data();
1203 namefull.Prepend("plots/pullsx_");
1204 namefull.Append(".eps");
1205 bx->Print(namefull.Data());
1206 namefull = name.Data();
1207 namefull.Prepend("plots/pullsx_");
1208 namefull.Append(".root");
1209 bx->Print(namefull.Data());
1211 namefull = name.Data();
1212 namefull.Prepend("plots/residualsz_");
1213 namefull.Append(".gif");
1214 az->Print(namefull.Data());
1215 namefull = name.Data();
1216 namefull.Prepend("plots/residualsz_");
1217 namefull.Append(".ps");
1218 az->Print(namefull.Data());
1219 namefull = name.Data();
1220 namefull.Prepend("plots/residualsz_");
1221 namefull.Append(".eps");
1222 az->Print(namefull.Data());
1223 namefull = name.Data();
1224 namefull.Prepend("plots/residualsz_");
1225 namefull.Append(".root");
1226 az->Print(namefull.Data());
1228 namefull = name.Data();
1229 namefull.Prepend("plots/pullsz_");
1230 namefull.Append(".gif");
1231 bz->Print(namefull.Data());
1232 namefull = name.Data();
1233 namefull.Prepend("plots/pullsz_");
1234 namefull.Append(".ps");
1235 bz->Print(namefull.Data());
1236 namefull = name.Data();
1237 namefull.Prepend("plots/pullsz_");
1238 namefull.Append(".eps");
1239 bz->Print(namefull.Data());
1240 namefull = name.Data();
1241 namefull.Prepend("plots/pullsz_");
1242 namefull.Append(".root");
1243 bz->Print(namefull.Data());
1245 namefull = name.Data();
1246 namefull.Prepend("plots/effVSntracklets_");
1247 namefull.Append(".gif");
1248 c1a->Print(namefull.Data());
1249 namefull = name.Data();
1250 namefull.Prepend("plots/effVSntracklets_");
1251 namefull.Append(".ps");
1252 c1a->Print(namefull.Data());
1253 namefull = name.Data();
1254 namefull.Prepend("plots/effVSntracklets_");
1255 namefull.Append(".eps");
1256 c1a->Print(namefull.Data());
1257 namefull = name.Data();
1258 namefull.Prepend("plots/effVSntracklets_");
1259 namefull.Append(".root");
1260 c1a->Print(namefull.Data());
1263 namefull = name.Data();
1264 namefull.Prepend("plots/effVSnitstracks_");
1265 namefull.Append(".gif");
1266 c1b->Print(namefull.Data());
1267 namefull = name.Data();
1268 namefull.Prepend("plots/effVSnitstracks_");
1269 namefull.Append(".ps");
1270 c1b->Print(namefull.Data());
1271 namefull = name.Data();
1272 namefull.Prepend("plots/effVSnitstracks_");
1273 namefull.Append(".eps");
1274 c1b->Print(namefull.Data());
1275 namefull = name.Data();
1276 namefull.Prepend("plots/effVSnitstracks_");
1277 namefull.Append(".root");
1278 c1b->Print(namefull.Data());
1281 namefull = name.Data();
1282 namefull.Prepend("plots/effVSnitstracks5or6_");
1283 namefull.Append(".gif");
1284 c1c->Print(namefull.Data());
1285 namefull = name.Data();
1286 namefull.Prepend("plots/effVSnitstracks5or6_");
1287 namefull.Append(".ps");
1288 c1c->Print(namefull.Data());
1289 namefull = name.Data();
1290 namefull.Prepend("plots/effVSnitstracks5or6_");
1291 namefull.Append(".eps");
1292 c1c->Print(namefull.Data());
1293 namefull = name.Data();
1294 namefull.Prepend("plots/effVSnitstracks5or6_");
1295 namefull.Append(".root");
1296 c1c->Print(namefull.Data());
1299 namefull = name.Data();
1300 namefull.Prepend("plots/effVSzMC_");
1301 namefull.Append(".gif");
1302 c1e->Print(namefull.Data());
1303 namefull = name.Data();
1304 namefull.Prepend("plots/effVSzMC_");
1305 namefull.Append(".ps");
1306 c1e->Print(namefull.Data());
1307 namefull = name.Data();
1308 namefull.Prepend("plots/effVSzMC_");
1309 namefull.Append(".eps");
1310 c1e->Print(namefull.Data());
1311 namefull = name.Data();
1312 namefull.Prepend("plots/effVSzMC_");
1313 namefull.Append(".root");
1314 c1e->Print(namefull.Data());
1318 namefull = name.Data();
1319 namefull.Prepend("plots/effVSntracklets2_");
1320 namefull.Append(".gif");
1321 c5->Print(namefull.Data());
1322 namefull = name.Data();
1323 namefull.Prepend("plots/effVSntracklets2_");
1324 namefull.Append(".ps");
1325 c5->Print(namefull.Data());
1326 namefull = name.Data();
1327 namefull.Prepend("plots/effVSntracklets2_");
1328 namefull.Append(".eps");
1329 c5->Print(namefull.Data());
1330 namefull = name.Data();
1331 namefull.Prepend("plots/effVSntracklets2_");
1332 namefull.Append(".root");
1333 c5->Print(namefull.Data());
1336 namefull = name.Data();
1337 namefull.Prepend("plots/sigmaVSntracklets_");
1338 namefull.Append(".gif");
1339 c6->Print(namefull.Data());
1340 namefull = name.Data();
1341 namefull.Prepend("plots/sigmaVSntracklets_");
1342 namefull.Append(".ps");
1343 c6->Print(namefull.Data());
1344 namefull = name.Data();
1345 namefull.Prepend("plots/sigmaVSntracklets_");
1346 namefull.Append(".eps");
1347 c6->Print(namefull.Data());
1348 namefull = name.Data();
1349 namefull.Prepend("plots/sigmaVSntracklets_");
1350 namefull.Append(".root");
1351 c6->Print(namefull.Data());
1354 namefull = name.Data();
1355 namefull.Prepend("plots/pullsVSntracklets_");
1356 namefull.Append(".gif");
1357 c8->Print(namefull.Data());
1358 namefull = name.Data();
1359 namefull.Prepend("plots/pullsVSntracklets_");
1360 namefull.Append(".ps");
1361 c8->Print(namefull.Data());
1362 namefull = name.Data();
1363 namefull.Prepend("plots/pullsVSntracklets_");
1364 namefull.Append(".eps");
1365 c8->Print(namefull.Data());
1366 namefull = name.Data();
1367 namefull.Prepend("plots/pullsVSntracklets_");
1368 namefull.Append(".root");
1369 c8->Print(namefull.Data());
1374 //----------------------------------------------------------------------------
1375 void TestRmTracks(Int_t iev=0) {
1378 delete gAlice->GetRunLoader();
1383 AliRunLoader *rl = AliRunLoader::Open("galice.root");
1385 cerr<<"Can not open session"<<endl;
1388 Int_t retval = rl->LoadgAlice();
1390 cerr<<"LoadgAlice returned error"<<endl;
1394 retval = rl->LoadHeader();
1396 cerr<<"LoadHeader returned error"<<endl;
1400 gAlice=rl->GetAliRun();
1401 // Get field from galice.root
1402 AliMagF *fiel = (AliMagF*)gAlice->Field();
1403 // Set the conversion constant between curvature and Pt
1404 AliTracker::SetFieldMap(fiel,kTRUE);
1406 Double_t truevtx[3];
1409 Double_t diffX,diffY,diffZ;
1410 Double_t diffXerr,diffYerr,diffZerr;
1411 Double_t multiplicity;
1414 // ----------- Create vertexer --------------------------
1415 AliESDVertex *initVertex = 0;
1416 TFile *invtx = new TFile("AliESDVertexMean.root");
1417 initVertex = (AliESDVertex*)invtx->Get("vtxmean");
1421 AliVertexerTracks *vertexer = new AliVertexerTracks();
1422 vertexer->SetVtxStart(initVertex);
1423 vertexer->SetOnlyFitter();
1424 vertexer->SetDebug(0); // set to 1 to see what it does
1426 Float_t diamondxy[2];
1427 diamondxy[0]=initVertex->GetXv();
1428 diamondxy[1]=initVertex->GetYv();
1429 //----------------------------------------------------------
1431 if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
1432 printf("AliESDs.root not found!\n");
1435 TFile *fin = TFile::Open("AliESDs.root");
1436 TTree *esdTree = (TTree*)fin->Get("esdTree");
1437 if(!esdTree) return;
1439 esdTree->SetBranchAddress("ESD",&esd);
1446 esdTree->GetEvent(iev);
1447 Int_t ntracks = esd->GetNumberOfTracks();
1448 printf("Number of tracks in ESD %d : %d\n",iev,ntracks);
1450 // example: do vertex without tracks 1 and 2 (AliESDtrack::GetID())
1452 // 1: tell vertexer to skip those two tracks
1458 vertexer->SetSkipTracks(nskip,skipped);
1459 AliESDVertex *vertex1 = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
1460 vertex1->PrintStatus();
1461 vertex1->PrintIndices();
1465 // 2: do vertex with all tracks
1469 vertexer->SetSkipTracks(nskip,skipped);
1470 AliESDVertex *vertex2 = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
1471 vertex2->PrintStatus();
1472 vertex2->PrintIndices();
1474 // 3: now remove those two tracks from vertex2
1476 TTree *trksTree = new TTree("trksTree","tree with tracks to be removed");
1477 AliESDtrack *esdTrack = 0;
1478 trksTree->Branch("tracks","AliESDtrack",&esdTrack);
1479 AliESDtrack *et = esd->GetTrack(1);
1480 esdTrack = new AliESDtrack(*et);
1483 AliESDtrack *et = esd->GetTrack(2);
1484 esdTrack = new AliESDtrack(*et);
1487 AliVertexerTracks vt;
1488 AliESDVertex *vertex3 = vt.RemoveTracksFromVertex(vertex2,trksTree,diamondxy);
1489 vertex3->PrintStatus();
1490 vertex3->PrintIndices();
1505 //----------------------------------------------------------------------------
1506 void AliComputeVtxMeanFromESD(TString file="AliESDs.root",
1507 Int_t mincontr=20) {
1508 //-----------------------------------------------------------------------
1509 // Compute weighted mean and cov. matrix from a set of AliESD
1510 //-----------------------------------------------------------------------
1512 gStyle->SetOptStat(0);
1513 //gStyle->SetOptFit(0);
1516 Double_t vtx[3],covvtx[6];
1519 AliESDVertex *vertex = 0;
1521 Int_t nc,events,total=0;
1530 Double_t sumwgtx=0.;
1531 Double_t sumwgty=0.;
1532 Double_t sumwgtz=0.;
1542 Double_t eavx,eavy,eavz,ewgtavx,ewgtavy,ewgtavz;
1544 TH1F* hx = new TH1F("hx","",200,-1000,1000);
1545 hx->SetXTitle("vertex x [#mu m]");
1546 hx->SetYTitle("events");
1547 TH1F* hy = new TH1F("hy","",200,-1000,1000);
1548 hy->SetXTitle("vertex y [#mu m]");
1549 hy->SetYTitle("events");
1550 TH1F* hz = new TH1F("hz","",200,-10,10);
1551 hz->SetXTitle("vertex z [cm]");
1552 hz->SetYTitle("events");
1555 for(Int_t i=1; i<100; i++) {
1556 TString inname="event.";
1559 inname.Append(file.Data());
1560 if(gSystem->AccessPathName(inname.Data())) continue;
1561 TFile *fin = TFile::Open(inname.Data());
1562 esdTree = (TTree*)fin->Get("esdTree");
1563 if(!esdTree) continue;
1564 esdTree->SetBranchAddress("ESD",&esd);
1565 events = esdTree->GetEntries();
1566 printf("Number of events in ESD tree: %d\n",events);
1569 for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
1570 esdTree->GetEvent(iev);
1571 vertex = (AliESDVertex*)esd->GetPrimaryVertex();
1572 if(!vertex) continue;
1573 nc = vertex->GetNContributors();
1574 if(nc<mincontr) continue;
1575 vertex->GetXYZ(vtx);
1576 vertex->GetSigmaXYZ(errvtx);
1577 //printf("%f %f %f %f %f %f\n",vtx[0],vtx[1],vtx[2],errvtx[0],errvtx[1],errvtx[2]);
1583 wgtavx += vtx[0]/errvtx[0]/errvtx[0];
1584 wgtavy += vtx[1]/errvtx[1]/errvtx[1];
1585 wgtavz += vtx[2]/errvtx[2]/errvtx[2];
1586 sumwgtx += 1./errvtx[0]/errvtx[0];
1587 sumwgty += 1./errvtx[1]/errvtx[1];
1588 sumwgtz += 1./errvtx[2]/errvtx[2];
1590 hx->Fill(1e4*vtx[0]);
1591 hy->Fill(1e4*vtx[1]);
1607 ewgtavx = 1./TMath::Sqrt(sumwgtx);
1608 ewgtavy = 1./TMath::Sqrt(sumwgty);
1609 ewgtavz = 1./TMath::Sqrt(sumwgtz);
1612 for(Int_t i=1; i<100; i++) {
1613 TString inname="event.";
1616 inname.Append(file.Data());
1617 if(gSystem->AccessPathName(inname.Data())) continue;
1618 TFile *fin = TFile::Open(inname.Data());
1619 esdTree = (TTree*)fin->Get("esdTree");
1620 if(!esdTree) continue;
1621 esdTree->SetBranchAddress("ESD",&esd);
1622 events = esdTree->GetEntries();
1623 for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
1624 esdTree->GetEvent(iev);
1625 vertex = (AliESDVertex*)esd->GetPrimaryVertex();
1626 if(!vertex) continue;
1627 nc = vertex->GetNContributors();
1628 if(nc<mincontr) continue;
1629 vertex->GetXYZ(vtx);
1630 vertex->GetSigmaXYZ(errvtx);
1632 varx += (vtx[0]-avx)*(vtx[0]-avx);
1633 vary += (vtx[1]-avy)*(vtx[1]-avy);
1634 varz += (vtx[2]-avz)*(vtx[2]-avz);
1635 covxy += (vtx[0]-avx)*(vtx[1]-avy);
1636 covxz += (vtx[0]-avx)*(vtx[2]-avz);
1637 covyz += (vtx[1]-avy)*(vtx[2]-avz);
1649 rmsx = TMath::Sqrt(varx);
1650 rmsy = TMath::Sqrt(vary);
1651 rmsz = TMath::Sqrt(varz);
1652 eavx = rmsx/TMath::Sqrt(sum);
1653 eavy = rmsy/TMath::Sqrt(sum);
1654 eavz = rmsz/TMath::Sqrt(sum);
1657 printf("\n\nNumber of events: Total %d, Used %d\n",total,sum);
1658 printf("Average:\n x = (%f +- %f) mum\n y = (%f +- %f) mum\n z = (%f +- %f) mu\n",1e4*avx,1e4*eavx,1e4*avy,1e4*eavy,1e4*avz,1e4*eavz);
1659 printf("Weighted Average:\n x = (%f +- %f) mum\n y = (%f +- %f) mum\n z = (%f +- %f) mum\n",1e4*wgtavx,1e4*ewgtavx,1e4*wgtavy,1e4*ewgtavy,1e4*wgtavz,1e4*ewgtavz);
1660 printf("RMS:\n x = %f mum\n y = %f mum\n z = %f mum\n",1e4*rmsx,1e4*rmsy,1e4*rmsz);
1662 TCanvas *c = new TCanvas("c","c",0,0,1000,500);
1666 TF1 *gx = new TF1("gx","gaus",-1000,1000);
1667 gx->SetLineColor(2);
1669 TF1 *gxx = (TF1*)gx->Clone("gxx");
1670 gxx->FixParameter(2,50.);
1671 gxx->SetLineStyle(2);
1675 TF1 *gy = new TF1("gy","gaus",-1000,1000);
1676 gy->SetLineColor(2);
1678 TF1 *gyy = (TF1*)gy->Clone("gyy");
1679 gyy->FixParameter(2,50.);
1680 gyy->SetLineStyle(2);
1684 TF1 *gz = new TF1("gz","gaus",-10,10);
1685 gz->SetLineColor(2);
1687 TF1 *gzz = (TF1*)gz->Clone("gzz");
1688 gzz->FixParameter(2,5.3);
1689 gzz->SetLineStyle(2);
1701 AliESDVertex *vtxmean = new AliESDVertex(vtx,covvtx,0.,(Int_t)sum);
1702 vtxmean->SetTitle("Vertex from weighted average");
1703 TFile *out = new TFile("AliESDVertexMean.root","recreate");
1704 vtxmean->Write("vtxmean");
1709 //----------------------------------------------------------------------------
1710 void AliComputeVtxMeanFromTree(TString file="AliESDVertices.root",
1711 Int_t mincontr=20) {
1712 //-----------------------------------------------------------------------
1713 // Compute weighted mean and cov. matrix from a set of AliESDVertices
1714 //-----------------------------------------------------------------------
1715 gStyle->SetOptStat(0);
1716 //gStyle->SetOptFit(0);
1719 Double_t vtx[3],covvtx[6];
1730 Double_t sumwgtx=0.;
1731 Double_t sumwgty=0.;
1732 Double_t sumwgtz=0.;
1736 Double_t eavx,eavy,eavz,ewgtavx,ewgtavy,ewgtavz;
1738 TH1F* hx = new TH1F("hx","",200,-1000,1000);
1739 hx->SetXTitle("vertex x [#mu m]");
1740 hx->SetYTitle("events");
1741 TH1F* hy = new TH1F("hy","",200,-1000,1000);
1742 hy->SetXTitle("vertex y [#mu m]");
1743 hy->SetYTitle("events");
1744 TH1F* hz = new TH1F("hz","",200,-10,10);
1745 hz->SetXTitle("vertex z [cm]");
1746 hz->SetYTitle("events");
1749 if(gSystem->AccessPathName(file.Data())) continue;
1750 TFile *fin = TFile::Open(file.Data());
1751 TTree *tree = (TTree*)fin->Get("TreeV");
1752 AliESDVertex *vertex = 0;
1753 tree->SetBranchAddress("AliESDVertex",&vertex);
1754 Int_t entries = (Int_t)tree->GetEntries();
1756 for(Int_t i=0; i<entries; i++) {
1759 nc = vertex->GetNContributors();
1760 if(nc<mincontr) continue;
1761 vertex->GetXYZ(vtx);
1762 vertex->GetSigmaXYZ(errvtx);
1763 //printf("%f %f %f %f %f %f\n",vtx[0],vtx[1],vtx[2],errvtx[0],errvtx[1],errvtx[2]);
1769 wgtavx += vtx[0]/errvtx[0]/errvtx[0];
1770 wgtavy += vtx[1]/errvtx[1]/errvtx[1];
1771 wgtavz += vtx[2]/errvtx[2]/errvtx[2];
1772 sumwgtx += 1./errvtx[0]/errvtx[0];
1773 sumwgty += 1./errvtx[1]/errvtx[1];
1774 sumwgtz += 1./errvtx[2]/errvtx[2];
1776 hx->Fill(1e4*vtx[0]);
1777 hy->Fill(1e4*vtx[1]);
1789 ewgtavx = 1./TMath::Sqrt(sumwgtx);
1790 ewgtavy = 1./TMath::Sqrt(sumwgty);
1791 ewgtavz = 1./TMath::Sqrt(sumwgtz);
1794 for(Int_t i=0; i<entries; i++) {
1796 nc = vertex->GetNContributors();
1797 if(nc<mincontr) continue;
1798 vertex->GetXYZ(vtx);
1799 vertex->GetSigmaXYZ(errvtx);
1801 varx += (vtx[0]-avx)*(vtx[0]-avx);
1802 vary += (vtx[1]-avy)*(vtx[1]-avy);
1803 varz += (vtx[2]-avz)*(vtx[2]-avz);
1804 covxy += (vtx[0]-avx)*(vtx[1]-avy);
1805 covxz += (vtx[0]-avx)*(vtx[2]-avz);
1806 covyz += (vtx[1]-avy)*(vtx[2]-avz);
1817 rmsx = TMath::Sqrt(varx);
1818 rmsy = TMath::Sqrt(vary);
1819 rmsz = TMath::Sqrt(varz);
1820 eavx = rmsx/TMath::Sqrt(sum);
1821 eavy = rmsy/TMath::Sqrt(sum);
1822 eavz = rmsz/TMath::Sqrt(sum);
1825 printf("\n\nNumber of events: Total %d, Used %d\n",total,sum);
1826 printf("Average:\n x = (%f +- %f) mum\n y = (%f +- %f) mum\n z = (%f +- %f) mu\n",1e4*avx,1e4*eavx,1e4*avy,1e4*eavy,1e4*avz,1e4*eavz);
1827 printf("Weighted Average:\n x = (%f +- %f) mum\n y = (%f +- %f) mum\n z = (%f +- %f) mum\n",1e4*wgtavx,1e4*ewgtavx,1e4*wgtavy,1e4*ewgtavy,1e4*wgtavz,1e4*ewgtavz);
1828 printf("RMS:\n x = %f mum\n y = %f mum\n z = %f mum\n",1e4*rmsx,1e4*rmsy,1e4*rmsz);
1830 TCanvas *c = new TCanvas("c","c",0,0,1000,500);
1834 TF1 *gx = new TF1("gx","gaus",-1000,1000);
1835 gx->SetLineColor(2);
1837 TF1 *gxx = (TF1*)gx->Clone("gxx");
1838 gxx->FixParameter(2,50.);
1839 gxx->SetLineStyle(2);
1843 TF1 *gy = new TF1("gy","gaus",-1000,1000);
1844 gy->SetLineColor(2);
1846 TF1 *gyy = (TF1*)gy->Clone("gyy");
1847 gyy->FixParameter(2,50.);
1848 gyy->SetLineStyle(2);
1852 TF1 *gz = new TF1("gz","gaus",-10,10);
1853 gz->SetLineColor(2);
1855 TF1 *gzz = (TF1*)gz->Clone("gzz");
1856 gzz->FixParameter(2,5.3);
1857 gzz->SetLineStyle(2);
1869 AliESDVertex *vtxmean = new AliESDVertex(vtx,covvtx,0.,(Int_t)sum);
1870 vtxmean->SetTitle("Vertex from weighted average");
1871 TFile *out = new TFile("AliESDVertexMean.root","recreate");
1872 vtxmean->Write("vtxmean");
1877 //---------------------------------------------------------------------------