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(AliTracker::GetBz());
91 vertexer->SetDebug(1); // set to 1 to see what it does
92 vertexer->SetVtxStart(initVertex);
93 if(!useMeanVtx) vertexer->SetConstraintOff();
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 AliESDEvent *esd = new AliESDEvent();
111 TTree *esdTree = (TTree*)fin->Get("esdTree");
113 esd->ReadFromTree(esdTree);
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 void VertexForOneEvent(Int_t iev=0,
255 Bool_t useMeanVtx=kFALSE) {
256 //------------------------------------------------------------------------
257 // Run vertex reconstruction for event iev
258 //------------------------------------------------------------------------
261 delete gAlice->GetRunLoader();
266 AliRunLoader *rl = AliRunLoader::Open("galice.root");
268 cerr<<"Can not open session"<<endl;
271 Int_t retval = rl->LoadgAlice();
273 cerr<<"LoadgAlice returned error"<<endl;
277 retval = rl->LoadHeader();
279 cerr<<"LoadHeader returned error"<<endl;
283 gAlice=rl->GetAliRun();
284 //rl->LoadKinematics();
285 // Get field from galice.root
286 AliMagF *fiel = (AliMagF*)gAlice->Field();
287 // Set the conversion constant between curvature and Pt
288 AliTracker::SetFieldMap(fiel,kTRUE);
293 Double_t diffX,diffY,diffZ;
294 Double_t diffXerr,diffYerr,diffZerr;
295 Double_t multiplicity;
298 // ----------- Create vertexer --------------------------
299 AliESDVertex *initVertex = 0;
301 // open the mean vertex
302 TFile *invtx = new TFile("AliESDVertexMean.root");
303 initVertex = (AliESDVertex*)invtx->Get("vtxmean");
307 Double_t pos[3]={0.,-0.,0.};
308 Double_t err[3]={3.,3.,5.};
309 initVertex = new AliESDVertex(pos,err);
311 AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
312 //vertexer->SetITSrefitNotRequired();
313 vertexer->SetVtxStart(initVertex);
314 if(!useMeanVtx) vertexer->SetConstraintOff();
315 //vertexer->SetMinTracks(2);
316 //vertexer->SetNSigmad0(nSigma);
317 //cout<<" Nsigma: "<<vertexer->GetNSigmad0()<<endl;
318 //vertexer->SetFinderAlgorithm(2);
319 //vertexer->SetDCAcut(0.1); // 1 mm
320 vertexer->SetDebug(1); // set to 1 to see what it does
321 //----------------------------------------------------------
323 if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
324 printf("AliESDs.root not found!\n");
327 TFile *fin = TFile::Open("AliESDs.root");
328 AliESDEvent *esd = new AliESDEvent();
329 TTree *esdTree = (TTree*)fin->Get("esdTree");
331 esd->ReadFromTree(esdTree);
337 //=================================================
338 // LOOK IN THE SIMULATION
339 // true vertex position
340 Int_t npart = (Int_t)gAlice->GetEvent(iev);
341 printf("particles %d\n",npart);
342 AliHeader *header = (AliHeader*)rl->GetHeader();
343 AliGenPythiaEventHeader *pyheader = (AliGenPythiaEventHeader*)header->GenEventHeader();
344 pyheader->PrimaryVertex(o);
348 printf("True Vertex: (%f, %f, %f)\n",o[0],o[1],o[2]);
350 esdTree->GetEvent(iev);
351 Int_t ntracks = esd->GetNumberOfTracks();
352 printf("Number of tracks in ESD %d : %d\n",iev,ntracks);
355 AliESDVertex *vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
356 nc = (Int_t)vertex->GetNContributors();
359 Int_t nitstracks = 0;
360 Int_t nitstracks1 = 0;
361 Int_t nitstracks2 = 0;
362 Int_t nitstracks3 = 0;
363 Int_t nitstracks4 = 0;
364 Int_t nitstracks5 = 0;
365 Int_t nitstracks6 = 0;
366 Int_t nitstracks5or6 = 0;
369 for(Int_t itrk=0; itrk<ntracks; itrk++) {
370 AliESDtrack *esdtrack = (AliESDtrack*)esd->GetTrack(itrk);
371 UInt_t status = esdtrack->GetStatus();
372 // only tracks found also in ITS
373 if(! (status&AliESDtrack::kITSrefit) ) continue;
375 Int_t npts = (Int_t)esdtrack->GetNcls(0);
376 if(npts==6) {nitstracks6++;nitstracks5or6++;}
377 if(npts==5) {nitstracks5++;nitstracks5or6++;}
378 if(npts==4) nitstracks4++;
379 if(npts==3) nitstracks3++;
380 if(npts==2) nitstracks2++;
381 if(npts==1) nitstracks1++;
383 printf("Number of kITSrefit tracks in ESD %d : %d\n",iev,nitstracks);
384 printf(" 6 pts: %d\n",nitstracks6);
385 printf(" 5 pts: %d\n",nitstracks5);
386 printf(" 4 pts: %d\n",nitstracks4);
387 printf(" 3 pts: %d\n",nitstracks3);
388 printf(" 2 pts: %d\n",nitstracks2);
389 printf(" 1 pts: %d\n",nitstracks1);
399 //----------------------------------------------------------------------------
400 Double_t FitFunc(Double_t *x,Double_t *par) {
401 return par[0]+par[1]/TMath::Sqrt(x[0]);
403 Int_t GetBin(Float_t mult) {
404 if(mult<5.) return 0;
405 if(mult<7.) return 1;
406 if(mult<12.) return 2;
407 if(mult<15.) return 3;
408 if(mult<22.) return 4;
411 Int_t GetBinTracklets(Float_t tracklets) {
412 if(tracklets<2.*5.) return 0;
413 if(tracklets<2.*7.) return 1;
414 if(tracklets<2.*12.) return 2;
415 if(tracklets<2.*15.) return 3;
416 if(tracklets<2.*22.) return 4;
419 Int_t GetBinZ(Float_t z) {
436 //----------------------------------------------------------------------------
437 void PlotVtxRes(TString inname="AliVertexerTracksTest.root") {
438 //---------------------------------------------------------------------
439 // Plot efficiency, resolutions, pulls
440 // [at least 10k pp events are needed]
441 //---------------------------------------------------------------------
445 Float_t ncminforzshift=1.;
447 gStyle->SetOptStat(0);
448 gStyle->SetOptFit(10001);
453 Float_t nESDTrks,nTrks,nTrks5or6,ntracklets;
454 Float_t diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,zMC,nc;
461 TH1F *hdx = new TH1F("hdx","",50,-600,600);
462 hdx->SetXTitle("#Delta x [#mu m]");
463 hdx->SetYTitle("events");
465 TH1F *hdy = new TH1F("hdy","",50,-600,600);
466 hdy->SetXTitle("#Delta y [#mu m]");
467 hdy->SetYTitle("events");
469 TH1F *hdz = new TH1F("hdz","",200,-600,600);
470 hdz->SetXTitle("#Delta z [#mu m]");
471 hdz->SetYTitle("events");
473 TH1F *hzMCVtx = new TH1F("hzMCVtx","events w/ vertex",30,-15,15);
474 hzMCVtx->SetXTitle("z_{true} [cm]");
475 hzMCVtx->SetYTitle("events");
478 TH1F *hzMCNoVtx = new TH1F("hzMCNoVtx","events w/o vertex",30,-15,15);
479 hzMCNoVtx->SetXTitle("z_{true} [cm]");
480 hzMCNoVtx->SetYTitle("events");
483 TH1F *hTrackletsVtx = new TH1F("hTrackletsVtx","events w/ vertex",51,-0.5,50.5);
484 hTrackletsVtx->SetXTitle("SPD tracklets");
485 hTrackletsVtx->SetYTitle("events");
487 TH1F *hTrackletsNoVtx = new TH1F("hTrackletsNoVtx","events w/o vertex",51,-0.5,50.5);
488 hTrackletsNoVtx->SetXTitle("SPD tracklets");
489 hTrackletsNoVtx->SetYTitle("events");
491 TH1F *hTracksVtx = new TH1F("hTracksVtx","events w/ vertex",51,-0.5,50.5);
492 hTracksVtx->SetXTitle("Number of reco tracks in TPC&ITS");
493 hTracksVtx->SetYTitle("events");
495 TH1F *hTracksNoVtx = new TH1F("hTracksNoVtx","events w/o vertex",51,-0.5,50.5);
496 hTracksNoVtx->SetXTitle("Number of reco tracks in TPC&ITS");
497 hTracksNoVtx->SetYTitle("events");
499 TH1F *hTracksVtx5or6 = new TH1F("hTracksVtx5or6","events w/ vertex",51,-0.5,50.5);
500 hTracksVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(cls>=5)");
501 hTracksVtx5or6->SetYTitle("events");
503 TH1F *hTracksNoVtx5or6 = new TH1F("hTracksNoVtx5or6","events w/o vertex",51,-0.5,50.5);
504 hTracksNoVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(cls>=5)");
505 hTracksNoVtx5or6->SetYTitle("events");
507 TH1F *hPullsx = new TH1F("hPullsx","",50,-7.,7.);
508 hPullsx->SetXTitle("#Delta x/#sqrt{<x,x>}");
509 hPullsx->SetYTitle("events");
511 TH1F *hPullsy = new TH1F("hPullsy","",50,-7.,7.);
512 hPullsy->SetXTitle("#Delta y/#sqrt{<y,y>}");
513 hPullsy->SetYTitle("events");
515 TH1F *hPullsz = new TH1F("hPullsz","",50,-7.,7.);
516 hPullsz->SetXTitle("#Delta z/#sqrt{<z,z>}");
517 hPullsz->SetYTitle("events");
519 TProfile *hntrackletsVSzMC = new TProfile("hntrackletsVSzMC","",30,-15,15,0,30);
520 hntrackletsVSzMC->SetXTitle("z_{true} [cm]");
521 hntrackletsVSzMC->SetYTitle("<SPD tracklets>");
523 TProfile *hnitstracksVSzMC = new TProfile("hnitstracksVSzMC","",30,-15,15,0,30);
524 hnitstracksVSzMC->SetXTitle("z_{true} [cm]");
525 hnitstracksVSzMC->SetYTitle("<tracks TPC&ITS>");
527 TProfile *hnitstracks5or6VSzMC = new TProfile("hnitstracks5or6VSzMC","",30,-15,15,0,30);
528 hnitstracks5or6VSzMC->SetXTitle("z_{true} [cm]");
529 hnitstracks5or6VSzMC->SetYTitle("<tracks TPC&ITS(cls>=5)>");
531 Double_t mult[6]={0,0,0,0,0,0};
532 Double_t mult2[6]={0,0,0,0,0,0};
533 Double_t sigmamult[6]={0,0,0,0,0,0};
534 Double_t sigmamult2[6]={0,0,0,0,0,0};
535 Double_t tottracks[6]={0,0,0,0,0,0};
536 Double_t fittrks[6]={0,0,0,0,0,0};
537 Double_t tracks[6]={0,0,0,0,0,0};
538 Double_t etracks[6]={0,0,0,0,0,0};
539 Double_t xres[6]={0,0,0,0,0,0};
540 Double_t exres[6]={0,0,0,0,0,0};
541 Double_t yres[6]={0,0,0,0,0,0};
542 Double_t eyres[6]={0,0,0,0,0,0};
543 Double_t zres[6]={0,0,0,0,0,0};
544 Double_t ezres[6]={0,0,0,0,0,0};
545 Double_t xpull[6]={0,0,0,0,0,0};
546 Double_t expull[6]={0,0,0,0,0,0};
547 Double_t ypull[6]={0,0,0,0,0,0};
548 Double_t eypull[6]={0,0,0,0,0,0};
549 Double_t zpull[6]={0,0,0,0,0,0};
550 Double_t ezpull[6]={0,0,0,0,0,0};
551 Int_t evts[6]={0,0,0,0,0,0};
552 Int_t totevts[6]={0,0,0,0,0,0};
553 Int_t vtx[6]={0,0,0,0,0,0};
554 Int_t all[6]={0,0,0,0,0,0};
555 Double_t probVtx[6]={0,0,0,0,0,0};
556 Double_t eprobVtx[6]={0,0,0,0,0,0};
561 TH1F *hdx0 = new TH1F("hdx0","x resolution - bin 0",50,-500,500);
562 TH1F *hdx1 = new TH1F("hdx1","x resolution - bin 1",50,-500,500);
563 TH1F *hdx2 = new TH1F("hdx2","x resolution - bin 2",50,-500,500);
564 TH1F *hdx3 = new TH1F("hdx3","x resolution - bin 3",50,-400,400);
565 TH1F *hdx4 = new TH1F("hdx4","x resolution - bin 4",50,-300,300);
566 TH1F *hdx5 = new TH1F("hdx5","x resolution - bin 5",50,-300,300);
567 TH1F *hdy0 = new TH1F("hdy0","y resolution - bin 0",50,-500,500);
568 TH1F *hdy1 = new TH1F("hdy1","y resolution - bin 1",50,-500,500);
569 TH1F *hdy2 = new TH1F("hdy2","y resolution - bin 2",50,-500,500);
570 TH1F *hdy3 = new TH1F("hdy3","y resolution - bin 3",50,-400,400);
571 TH1F *hdy4 = new TH1F("hdy4","y resolution - bin 4",50,-300,300);
572 TH1F *hdy5 = new TH1F("hdy5","y resolution - bin 5",50,-300,300);
573 TH1F *hdz0 = new TH1F("hdz0","z resolution - bin 0",50,-1000,1000);
574 TH1F *hdz1 = new TH1F("hdz1","z resolution - bin 1",50,-800,800);
575 TH1F *hdz2 = new TH1F("hdz2","z resolution - bin 2",50,-800,800);
576 TH1F *hdz3 = new TH1F("hdz3","z resolution - bin 3",50,-600,600);
577 TH1F *hdz4 = new TH1F("hdz4","z resolution - bin 4",50,-500,500);
578 TH1F *hdz5 = new TH1F("hdz5","z resolution - bin 5",50,-500,500);
580 TH1F *hdz_z = NULL; hdz_z = new TH1F[15];
581 for(Int_t l=0;l<15;l++) hdz_z[l] = *hdz;
584 TH1F *hpx0 = new TH1F("hpx0","x pull - bin 0",50,-7,7);
585 TH1F *hpx1 = new TH1F("hpx1","x pull - bin 1",50,-7,7);
586 TH1F *hpx2 = new TH1F("hpx2","x pull - bin 2",50,-7,7);
587 TH1F *hpx3 = new TH1F("hpx3","x pull - bin 3",50,-7,7);
588 TH1F *hpx4 = new TH1F("hpx4","x pull - bin 4",50,-7,7);
589 TH1F *hpx5 = new TH1F("hpx5","x pull - bin 5",50,-7,7);
590 TH1F *hpy0 = new TH1F("hpy0","y pull - bin 0",50,-7,7);
591 TH1F *hpy1 = new TH1F("hpy1","y pull - bin 1",50,-7,7);
592 TH1F *hpy2 = new TH1F("hpy2","y pull - bin 2",50,-7,7);
593 TH1F *hpy3 = new TH1F("hpy3","y pull - bin 3",50,-7,7);
594 TH1F *hpy4 = new TH1F("hpy4","y pull - bin 4",50,-7,7);
595 TH1F *hpy5 = new TH1F("hpy5","y pull - bin 5",50,-7,7);
596 TH1F *hpz0 = new TH1F("hpz0","z pull - bin 0",50,-7,7);
597 TH1F *hpz1 = new TH1F("hpz1","z pull - bin 1",50,-7,7);
598 TH1F *hpz2 = new TH1F("hpz2","z pull - bin 2",50,-7,7);
599 TH1F *hpz3 = new TH1F("hpz3","z pull - bin 3",50,-7,7);
600 TH1F *hpz4 = new TH1F("hpz4","z pull - bin 4",50,-7,7);
601 TH1F *hpz5 = new TH1F("hpz5","z pull - bin 5",50,-7,7);
604 TH1F *hmult = new TH1F("hmult","hmult",50,-0.5,49.5);
606 TFile *in = new TFile(inname.Data());
607 TNtuple *nt = (TNtuple*)in->Get("ntVtxRes");
608 nt->SetBranchAddress("triggered",&triggered);
609 nt->SetBranchAddress("ntracks",&nESDTrks);
610 nt->SetBranchAddress("nitstracks",&nTrks);
611 nt->SetBranchAddress("nitstracks5or6",&nTrks5or6);
612 nt->SetBranchAddress("diffX",&diffX);
613 nt->SetBranchAddress("diffY",&diffY);
614 nt->SetBranchAddress("diffZ",&diffZ);
615 nt->SetBranchAddress("diffXerr",&diffXerr);
616 nt->SetBranchAddress("diffYerr",&diffYerr);
617 nt->SetBranchAddress("diffZerr",&diffZerr);
618 nt->SetBranchAddress("multiplicity",&ntracklets);
619 nt->SetBranchAddress("zMC",&zMC);
620 nt->SetBranchAddress("nc",&nc);
621 Int_t entries = (Int_t)nt->GetEntries();
624 for(Int_t i=0; i<entries; i++) {
625 nbytes += nt->GetEvent(i);
626 // only triggered events
627 if(triggered<0.5) continue;
628 // consider only events with zMCmin<|zMC|<zMCmax
629 if(TMath::Abs(zMC)<zMCmin || TMath::Abs(zMC)>zMCmax) continue;
631 bin = GetBinTracklets(ntracklets);
633 mult[bin] += ntracklets;
634 hmult->Fill(ntracklets);
636 hntrackletsVSzMC->Fill(zMC,ntracklets);
637 hnitstracksVSzMC->Fill(zMC,nTrks);
638 hnitstracks5or6VSzMC->Fill(zMC,nTrks5or6);
639 if(TMath::Abs(diffX) < 1000.) { // vtx found - closer than 1 mm
641 mult2[bin] += ntracklets;
643 tottracks[bin] += nTrks;
645 if(bin==0) hdx0->Fill(diffX);
646 if(bin==1) hdx1->Fill(diffX);
647 if(bin==2) hdx2->Fill(diffX);
648 if(bin==3) hdx3->Fill(diffX);
649 if(bin==4) hdx4->Fill(diffX);
650 if(bin==5) hdx5->Fill(diffX);
651 if(bin==0) hdy0->Fill(diffY);
652 if(bin==1) hdy1->Fill(diffY);
653 if(bin==2) hdy2->Fill(diffY);
654 if(bin==3) hdy3->Fill(diffY);
655 if(bin==4) hdy4->Fill(diffY);
656 if(bin==5) hdy5->Fill(diffY);
657 if(bin==0) hdz0->Fill(diffZ);
658 if(bin==1) hdz1->Fill(diffZ);
659 if(bin==2) hdz2->Fill(diffZ);
660 if(bin==3) hdz3->Fill(diffZ);
661 if(bin==4) hdz4->Fill(diffZ);
662 if(bin==5) hdz5->Fill(diffZ);
666 if(nc>=ncminforzshift) hdz_z[zbin].Fill(diffZ);
668 hTrackletsVtx->Fill(ntracklets);
669 hTracksVtx->Fill(nTrks);
670 hTracksVtx5or6->Fill(nTrks5or6);
671 hPullsx->Fill(diffXerr);
672 hPullsy->Fill(diffYerr);
673 hPullsz->Fill(diffZerr);
674 if(bin==0) hpx0->Fill(diffXerr);
675 if(bin==1) hpx1->Fill(diffXerr);
676 if(bin==2) hpx2->Fill(diffXerr);
677 if(bin==3) hpx3->Fill(diffXerr);
678 if(bin==4) hpx4->Fill(diffXerr);
679 if(bin==5) hpx5->Fill(diffXerr);
680 if(bin==0) hpy0->Fill(diffYerr);
681 if(bin==1) hpy1->Fill(diffYerr);
682 if(bin==2) hpy2->Fill(diffYerr);
683 if(bin==3) hpy3->Fill(diffYerr);
684 if(bin==4) hpy4->Fill(diffYerr);
685 if(bin==5) hpy5->Fill(diffYerr);
686 if(bin==0) hpz0->Fill(diffZerr);
687 if(bin==1) hpz1->Fill(diffZerr);
688 if(bin==2) hpz2->Fill(diffZerr);
689 if(bin==3) hpz3->Fill(diffZerr);
690 if(bin==4) hpz4->Fill(diffZerr);
691 if(bin==5) hpz5->Fill(diffZerr);
694 hzMCNoVtx->Fill(zMC);
695 hTrackletsNoVtx->Fill(ntracklets);
696 hTracksNoVtx->Fill(nTrks);
697 hTracksNoVtx5or6->Fill(nTrks5or6);
703 for(Int_t k=0;k<6;k++) {
705 mult[k] /= totevts[k];
706 tottracks[k] /= evts[k];
707 probVtx[k] = (Double_t)vtx[k]/all[k];
708 eprobVtx[k] = (Double_t)vtx[k]/all[k]/all[k]+(Double_t)vtx[k]*vtx[k]/all[k]/all[k]/all[k];
709 eprobVtx[k] = TMath::Sqrt(eprobVtx[k]);
712 // calculate spread in multiplicity
713 for(Int_t i=0; i<entries; i++) {
714 nbytes += nt->GetEvent(i);
715 bin = GetBinTracklets(ntracklets);
716 sigmamult[bin] += (ntracklets-mult[bin])*(ntracklets-mult[bin])/totevts[bin];
717 if(diffX < 90000.) sigmamult2[bin] += (ntracklets-mult2[bin])*(ntracklets-mult2[bin])/evts[bin];
720 for(Int_t k=0;k<6;k++) {
721 sigmamult[k] = TMath::Sqrt(sigmamult[k]);
722 sigmamult2[k] = TMath::Sqrt(sigmamult2[k]);
725 // fraction of found vertices
726 Float_t frac = (Float_t)nEvVtx/(nEvVtx+nEvNoVtx);
727 printf(" Fraction of events with vertex: %f\n",frac);
731 TF1 *g = new TF1("g","gaus");
733 TF1 *line = new TF1("line","pol1");
734 line->SetLineWidth(3);
735 TF1 *func = new TF1("func",FitFunc,2,80,2);
736 func->SetLineWidth(1);
738 Double_t x1,y1,sigmafit;
744 gStyle->SetOptFit(1111);
747 TCanvas *c0 = new TCanvas("c0","c0",0,0,1000,400);
750 hntrackletsVSzMC->SetMinimum(0);
751 hntrackletsVSzMC->SetMaximum(14);
752 hntrackletsVSzMC->Draw("profs");
754 hnitstracksVSzMC->SetMinimum(0);
755 hnitstracksVSzMC->SetMaximum(14);
756 hnitstracksVSzMC->Draw("profs");
758 hnitstracks5or6VSzMC->SetMinimum(0);
759 hnitstracks5or6VSzMC->SetMaximum(14);
760 hnitstracks5or6VSzMC->Draw("profs");
763 // Tracks in ITS for events w/ and w/o vertex
764 TCanvas *c1 = new TCanvas("c1","c1",0,0,1000,500);
769 hTracksNoVtx->Draw();
771 // probability vs ntracklets
772 TCanvas *c1a = new TCanvas("c1a","c1a",0,0,500,500);
773 TH1F *hTotTracklets = (TH1F*)hTrackletsVtx->Clone("hTotTracklets");
774 hTotTracklets->Add(hTrackletsNoVtx);
775 TH1F *hProbTracklets = (TH1F*)hTrackletsVtx->Clone("hProbTracklets");
776 hProbTracklets->Divide(hTotTracklets);
777 hProbTracklets->SetYTitle("Probability");
778 hProbTracklets->SetTitle("Probability to find the vertex");
779 hProbTracklets->Draw();
781 // probability vs ntracks
782 TCanvas *c1b = new TCanvas("c1b","c1b",0,0,500,500);
783 TH1F *hTotTracks = (TH1F*)hTracksVtx->Clone("hTotTracks");
784 hTotTracks->Add(hTracksNoVtx);
785 TH1F *hProbTracks = (TH1F*)hTracksVtx->Clone("hProbTracks");
786 hProbTracks->Divide(hTotTracks);
787 hProbTracks->SetYTitle("Probability");
788 hProbTracks->SetTitle("Probability to find the vertex");
791 TCanvas *c1c = new TCanvas("c1c","c1c",0,0,500,500);
792 TH1F *hTotTracks5or6 = (TH1F*)hTracksVtx5or6->Clone("hTotTracks5or6");
793 hTotTracks5or6->Add(hTracksNoVtx5or6);
794 TH1F *hProbTracks5or6 = (TH1F*)hTracksVtx5or6->Clone("hProbTracks5or6");
795 hProbTracks5or6->Divide(hTotTracks5or6);
796 hProbTracks5or6->SetYTitle("Probability");
797 hProbTracks5or6->SetTitle("Probability to find the vertex");
798 hProbTracks5or6->Draw();
801 // probability vs zMC
802 TCanvas *c1e = new TCanvas("c1e","c1e",0,0,500,500);
803 TH1F *hTotzMC = (TH1F*)hzMCVtx->Clone("hTotzMC");
804 hTotzMC->Add(hzMCNoVtx);
805 TH1F *hProbzMC = (TH1F*)hzMCVtx->Clone("hProbzMC");
806 hProbzMC->Divide(hTotzMC);
807 hProbzMC->SetYTitle("Probability");
808 hProbzMC->SetTitle("Probability to find the vertex");
811 // Global resolutions
812 TCanvas *c2 = new TCanvas("c2","c2",0,0,1000,400);
816 g->SetRange(-1.*hdx->GetRMS(),+1.*hdx->GetRMS());
820 g->SetRange(-1.*hdy->GetRMS(),+1.*hdy->GetRMS());
824 g->SetRange(-1.*hdz->GetRMS(),+1.*hdz->GetRMS());
828 TCanvas *c4 = new TCanvas("c4","c4",0,0,1000,400);
832 g->SetRange(-3.*hPullsx->GetRMS(),+3.*hPullsx->GetRMS());
833 hPullsx->Fit("g","R");
836 g->SetRange(-3.*hPullsy->GetRMS(),+3.*hPullsy->GetRMS());
837 hPullsy->Fit("g","R");
840 g->SetRange(-3.*hPullsz->GetRMS(),+3.*hPullsz->GetRMS());
841 hPullsz->Fit("g","R");
843 // probability VS multiplicity
844 TCanvas *c5 = new TCanvas("c5","c5",0,0,600,500);
845 TH2F *fr5 = new TH2F("fr5","",2,0,80,2,0,1.1);
846 fr5->SetXTitle("SPD tracklets");
847 fr5->SetYTitle("Probability");
849 TGraphErrors *gr5 = new TGraphErrors(6,mult,probVtx,sigmamult,eprobVtx);
851 gr5->SetMarkerStyle(21);
854 // resolutions VS multiplicity
857 TCanvas *ax = new TCanvas("ax","ax",0,0,900,700);
861 g->SetRange(-3.*hdx0->GetRMS(),+3.*hdx0->GetRMS());
863 xres[0]=g->GetParameter(2);
864 exres[0]=g->GetParError(2);
867 g->SetRange(-3.*hdx1->GetRMS(),+3.*hdx1->GetRMS());
869 xres[1]=g->GetParameter(2);
870 exres[1]=g->GetParError(2);
873 g->SetRange(-3.*hdx2->GetRMS(),+3.*hdx2->GetRMS());
875 xres[2]=g->GetParameter(2);
876 exres[2]=g->GetParError(2);
879 g->SetRange(-3.*hdx3->GetRMS(),+3.*hdx3->GetRMS());
881 xres[3]=g->GetParameter(2);
882 exres[3]=g->GetParError(2);
885 g->SetRange(-3.*hdx4->GetRMS(),+3.*hdx4->GetRMS());
887 xres[4]=g->GetParameter(2);
888 exres[4]=g->GetParError(2);
891 g->SetRange(-3.*hdx5->GetRMS(),+3.*hdx5->GetRMS());
893 xres[5]=g->GetParameter(2);
894 exres[5]=g->GetParError(2);
896 TCanvas *ay = new TCanvas("ay","ay",0,0,900,700);
900 g->SetRange(-3.*hdy0->GetRMS(),+3.*hdy0->GetRMS());
902 yres[0]=g->GetParameter(2);
903 eyres[0]=g->GetParError(2);
906 g->SetRange(-3.*hdy1->GetRMS(),+3.*hdy1->GetRMS());
908 yres[1]=g->GetParameter(2);
909 eyres[1]=g->GetParError(2);
912 g->SetRange(-3.*hdy2->GetRMS(),+3.*hdy2->GetRMS());
914 yres[2]=g->GetParameter(2);
915 eyres[2]=g->GetParError(2);
918 g->SetRange(-3.*hdy3->GetRMS(),+3.*hdy3->GetRMS());
920 yres[3]=g->GetParameter(2);
921 eyres[3]=g->GetParError(2);
924 g->SetRange(-3.*hdy4->GetRMS(),+3.*hdy4->GetRMS());
926 yres[4]=g->GetParameter(2);
927 eyres[4]=g->GetParError(2);
930 g->SetRange(-3.*hdy5->GetRMS(),+3.*hdy5->GetRMS());
932 yres[5]=g->GetParameter(2);
933 eyres[5]=g->GetParError(2);
935 TCanvas *az = new TCanvas("az","az",0,0,900,700);
939 g->SetRange(-3.*hdz0->GetRMS(),+3.*hdz0->GetRMS());
941 zres[0]=g->GetParameter(2);
942 ezres[0]=g->GetParError(2);
945 g->SetRange(-3.*hdz1->GetRMS(),+3.*hdz1->GetRMS());
947 zres[1]=g->GetParameter(2);
948 ezres[1]=g->GetParError(2);
951 g->SetRange(-3.*hdz2->GetRMS(),+3.*hdz2->GetRMS());
953 zres[2]=g->GetParameter(2);
954 ezres[2]=g->GetParError(2);
957 g->SetRange(-3.*hdz3->GetRMS(),+3.*hdz3->GetRMS());
959 zres[3]=g->GetParameter(2);
960 ezres[3]=g->GetParError(2);
963 g->SetRange(-3.*hdz4->GetRMS(),+3.*hdz4->GetRMS());
965 zres[4]=g->GetParameter(2);
966 ezres[4]=g->GetParError(2);
969 g->SetRange(-3.*hdz5->GetRMS(),+3.*hdz5->GetRMS());
971 zres[5]=g->GetParameter(2);
972 ezres[5]=g->GetParError(2);
976 // pulls VS multiplicity
979 TCanvas *bx = new TCanvas("bx","bx",0,0,900,700);
983 g->SetRange(-3.*hpx0->GetRMS(),+3.*hpx0->GetRMS());
985 xpull[0]=g->GetParameter(2);
986 expull[0]=g->GetParError(2);
989 g->SetRange(-3.*hpx1->GetRMS(),+3.*hpx1->GetRMS());
991 xpull[1]=g->GetParameter(2);
992 expull[1]=g->GetParError(2);
995 g->SetRange(-3.*hpx2->GetRMS(),+3.*hpx2->GetRMS());
997 xpull[2]=g->GetParameter(2);
998 expull[2]=g->GetParError(2);
1001 g->SetRange(-3.*hpx3->GetRMS(),+3.*hpx3->GetRMS());
1003 xpull[3]=g->GetParameter(2);
1004 expull[3]=g->GetParError(2);
1007 g->SetRange(-3.*hpx4->GetRMS(),+3.*hpx4->GetRMS());
1009 xpull[4]=g->GetParameter(2);
1010 expull[4]=g->GetParError(2);
1013 g->SetRange(-3.*hpx5->GetRMS(),+3.*hpx5->GetRMS());
1015 xpull[5]=g->GetParameter(2);
1016 expull[5]=g->GetParError(2);
1018 TCanvas *by = new TCanvas("by","by",0,0,900,700);
1022 g->SetRange(-3.*hpy0->GetRMS(),+3.*hpy0->GetRMS());
1024 ypull[0]=g->GetParameter(2);
1025 eypull[0]=g->GetParError(2);
1028 g->SetRange(-3.*hpy1->GetRMS(),+3.*hpy1->GetRMS());
1030 ypull[1]=g->GetParameter(2);
1031 eypull[1]=g->GetParError(2);
1034 g->SetRange(-3.*hpy2->GetRMS(),+3.*hpy2->GetRMS());
1036 ypull[2]=g->GetParameter(2);
1037 eypull[2]=g->GetParError(2);
1040 g->SetRange(-3.*hpy3->GetRMS(),+3.*hpy3->GetRMS());
1042 ypull[3]=g->GetParameter(2);
1043 eypull[3]=g->GetParError(2);
1046 g->SetRange(-3.*hpy4->GetRMS(),+3.*hpy4->GetRMS());
1048 ypull[4]=g->GetParameter(2);
1049 eypull[4]=g->GetParError(2);
1052 g->SetRange(-3.*hpy5->GetRMS(),+3.*hpy5->GetRMS());
1054 ypull[5]=g->GetParameter(2);
1055 eypull[5]=g->GetParError(2);
1057 TCanvas *bz = new TCanvas("bz","bz",0,0,900,700);
1061 g->SetRange(-3.*hpz0->GetRMS(),+3.*hpz0->GetRMS());
1063 zpull[0]=g->GetParameter(2);
1064 ezpull[0]=g->GetParError(2);
1067 g->SetRange(-3.*hpz1->GetRMS(),+3.*hpz1->GetRMS());
1069 zpull[1]=g->GetParameter(2);
1070 ezpull[1]=g->GetParError(2);
1073 g->SetRange(-3.*hpz2->GetRMS(),+3.*hpz2->GetRMS());
1075 zpull[2]=g->GetParameter(2);
1076 ezpull[2]=g->GetParError(2);
1079 g->SetRange(-3.*hpz3->GetRMS(),+3.*hpz3->GetRMS());
1081 zpull[3]=g->GetParameter(2);
1082 ezpull[3]=g->GetParError(2);
1085 g->SetRange(-3.*hpz4->GetRMS(),+3.*hpz4->GetRMS());
1087 zpull[4]=g->GetParameter(2);
1088 ezpull[4]=g->GetParError(2);
1091 g->SetRange(-3.*hpz5->GetRMS(),+3.*hpz5->GetRMS());
1093 zpull[5]=g->GetParameter(2);
1094 ezpull[5]=g->GetParError(2);
1096 gStyle->SetOptFit(0);
1098 // res VS ntracklets
1099 TCanvas *c6 = new TCanvas("c6","c6",0,0,600,500);
1100 TH2F *fr6 = new TH2F("fr6","",2,0,80,2,0,240);
1101 fr6->SetXTitle("SPD tracklets");
1102 fr6->SetYTitle("#sigma [#mu m]");
1104 sigmamult2[0]=sigmamult2[1]=sigmamult2[2]=sigmamult2[3]=sigmamult2[4]=sigmamult2[5]=0.;
1105 TGraphErrors *gr6x = new TGraphErrors(6,mult2,xres,sigmamult2,exres);
1107 gr6x->SetMarkerStyle(22);
1108 //gr6x->Fit("func","E,R");
1109 TGraphErrors *gr6z = new TGraphErrors(6,mult2,zres,sigmamult2,ezres);
1111 gr6z->SetMarkerStyle(26);
1112 //gr6z->Fit("func","E,R");
1113 TLegend *leg6 = new TLegend(0.6,0.75,0.9,0.9);
1114 leg6->AddEntry(gr6x,"x/y coordinate","p");
1115 leg6->AddEntry(gr6z,"z coordinate","p");
1116 leg6->SetFillStyle(0);
1117 leg6->SetBorderSize(0);
1120 // pull VS ntracklets
1121 TCanvas *c8 = new TCanvas("c8","c8",0,0,600,500);
1122 TH2F *fr8 = new TH2F("fr8","",2,0,80,2,0,2.);
1123 fr8->SetXTitle("SPD tracklets");
1124 fr8->SetYTitle("pull");
1126 TGraphErrors *gr8x = new TGraphErrors(6,mult2,xpull,sigmamult2,expull);
1128 gr8x->SetMarkerStyle(22);
1129 TGraphErrors *gr8z = new TGraphErrors(6,mult2,zpull,sigmamult2,ezpull);
1131 gr8z->SetMarkerStyle(26);
1132 TLegend *leg8 = new TLegend(0.6,0.75,0.9,0.9);
1133 leg8->AddEntry(gr6x,"x/y coordinate","p");
1134 leg8->AddEntry(gr6z,"z coordinate","p");
1135 leg8->SetFillStyle(0);
1136 leg8->SetBorderSize(0);
1138 TLine *l8 = new TLine(0,1,80,1);
1139 l8->SetLineStyle(2);
1143 // mean z res VS zMC
1144 Float_t zmc[15]={-14,-12,-10,-8,-6,-4,-2,0,2,4,6,8,10,12,14};
1146 Float_t dzmean[15],edzmean[15];
1148 TCanvas *zz = new TCanvas("zz","zz",0,0,1000,1000);
1153 g->SetRange(-3.*hdz_z[l].GetRMS(),+3.*hdz_z[l].GetRMS());
1154 hdz_z[l].Fit("g","Q");
1155 dzmean[l]=g->GetParameter(1);
1156 edzmean[l]=g->GetParError(1);
1159 TCanvas *zzz = new TCanvas("zzz","zzz",0,0,600,500);
1160 TH2F *fr9 = new TH2F("fr9","",2,-15,15,2,-50,50);
1161 fr9->SetXTitle("z_{true} [cm]");
1162 fr9->SetYTitle("<z_{rec} - z_{true}> [#mu m]");
1164 TGraphErrors *grz = new TGraphErrors(15,zmc,dzmean,ezmc,edzmean);
1166 grz->SetMarkerStyle(20);
1171 //----------------------------------------------------------------------------
1172 void SaveFigures(TString name="dummy") {
1176 namefull = name.Data();
1177 namefull.Prepend("plots/residualsx_");
1178 namefull.Append(".gif");
1179 ax->Print(namefull.Data());
1180 namefull = name.Data();
1181 namefull.Prepend("plots/residualsx_");
1182 namefull.Append(".ps");
1183 ax->Print(namefull.Data());
1184 namefull = name.Data();
1185 namefull.Prepend("plots/residualsx_");
1186 namefull.Append(".eps");
1187 ax->Print(namefull.Data());
1188 namefull = name.Data();
1189 namefull.Prepend("plots/residualsx_");
1190 namefull.Append(".root");
1191 ax->Print(namefull.Data());
1193 namefull = name.Data();
1194 namefull.Prepend("plots/pullsx_");
1195 namefull.Append(".gif");
1196 bx->Print(namefull.Data());
1197 namefull = name.Data();
1198 namefull.Prepend("plots/pullsx_");
1199 namefull.Append(".ps");
1200 bx->Print(namefull.Data());
1201 namefull = name.Data();
1202 namefull.Prepend("plots/pullsx_");
1203 namefull.Append(".eps");
1204 bx->Print(namefull.Data());
1205 namefull = name.Data();
1206 namefull.Prepend("plots/pullsx_");
1207 namefull.Append(".root");
1208 bx->Print(namefull.Data());
1210 namefull = name.Data();
1211 namefull.Prepend("plots/residualsz_");
1212 namefull.Append(".gif");
1213 az->Print(namefull.Data());
1214 namefull = name.Data();
1215 namefull.Prepend("plots/residualsz_");
1216 namefull.Append(".ps");
1217 az->Print(namefull.Data());
1218 namefull = name.Data();
1219 namefull.Prepend("plots/residualsz_");
1220 namefull.Append(".eps");
1221 az->Print(namefull.Data());
1222 namefull = name.Data();
1223 namefull.Prepend("plots/residualsz_");
1224 namefull.Append(".root");
1225 az->Print(namefull.Data());
1227 namefull = name.Data();
1228 namefull.Prepend("plots/pullsz_");
1229 namefull.Append(".gif");
1230 bz->Print(namefull.Data());
1231 namefull = name.Data();
1232 namefull.Prepend("plots/pullsz_");
1233 namefull.Append(".ps");
1234 bz->Print(namefull.Data());
1235 namefull = name.Data();
1236 namefull.Prepend("plots/pullsz_");
1237 namefull.Append(".eps");
1238 bz->Print(namefull.Data());
1239 namefull = name.Data();
1240 namefull.Prepend("plots/pullsz_");
1241 namefull.Append(".root");
1242 bz->Print(namefull.Data());
1244 namefull = name.Data();
1245 namefull.Prepend("plots/effVSntracklets_");
1246 namefull.Append(".gif");
1247 c1a->Print(namefull.Data());
1248 namefull = name.Data();
1249 namefull.Prepend("plots/effVSntracklets_");
1250 namefull.Append(".ps");
1251 c1a->Print(namefull.Data());
1252 namefull = name.Data();
1253 namefull.Prepend("plots/effVSntracklets_");
1254 namefull.Append(".eps");
1255 c1a->Print(namefull.Data());
1256 namefull = name.Data();
1257 namefull.Prepend("plots/effVSntracklets_");
1258 namefull.Append(".root");
1259 c1a->Print(namefull.Data());
1262 namefull = name.Data();
1263 namefull.Prepend("plots/effVSnitstracks_");
1264 namefull.Append(".gif");
1265 c1b->Print(namefull.Data());
1266 namefull = name.Data();
1267 namefull.Prepend("plots/effVSnitstracks_");
1268 namefull.Append(".ps");
1269 c1b->Print(namefull.Data());
1270 namefull = name.Data();
1271 namefull.Prepend("plots/effVSnitstracks_");
1272 namefull.Append(".eps");
1273 c1b->Print(namefull.Data());
1274 namefull = name.Data();
1275 namefull.Prepend("plots/effVSnitstracks_");
1276 namefull.Append(".root");
1277 c1b->Print(namefull.Data());
1280 namefull = name.Data();
1281 namefull.Prepend("plots/effVSnitstracks5or6_");
1282 namefull.Append(".gif");
1283 c1c->Print(namefull.Data());
1284 namefull = name.Data();
1285 namefull.Prepend("plots/effVSnitstracks5or6_");
1286 namefull.Append(".ps");
1287 c1c->Print(namefull.Data());
1288 namefull = name.Data();
1289 namefull.Prepend("plots/effVSnitstracks5or6_");
1290 namefull.Append(".eps");
1291 c1c->Print(namefull.Data());
1292 namefull = name.Data();
1293 namefull.Prepend("plots/effVSnitstracks5or6_");
1294 namefull.Append(".root");
1295 c1c->Print(namefull.Data());
1298 namefull = name.Data();
1299 namefull.Prepend("plots/effVSzMC_");
1300 namefull.Append(".gif");
1301 c1e->Print(namefull.Data());
1302 namefull = name.Data();
1303 namefull.Prepend("plots/effVSzMC_");
1304 namefull.Append(".ps");
1305 c1e->Print(namefull.Data());
1306 namefull = name.Data();
1307 namefull.Prepend("plots/effVSzMC_");
1308 namefull.Append(".eps");
1309 c1e->Print(namefull.Data());
1310 namefull = name.Data();
1311 namefull.Prepend("plots/effVSzMC_");
1312 namefull.Append(".root");
1313 c1e->Print(namefull.Data());
1317 namefull = name.Data();
1318 namefull.Prepend("plots/effVSntracklets2_");
1319 namefull.Append(".gif");
1320 c5->Print(namefull.Data());
1321 namefull = name.Data();
1322 namefull.Prepend("plots/effVSntracklets2_");
1323 namefull.Append(".ps");
1324 c5->Print(namefull.Data());
1325 namefull = name.Data();
1326 namefull.Prepend("plots/effVSntracklets2_");
1327 namefull.Append(".eps");
1328 c5->Print(namefull.Data());
1329 namefull = name.Data();
1330 namefull.Prepend("plots/effVSntracklets2_");
1331 namefull.Append(".root");
1332 c5->Print(namefull.Data());
1335 namefull = name.Data();
1336 namefull.Prepend("plots/sigmaVSntracklets_");
1337 namefull.Append(".gif");
1338 c6->Print(namefull.Data());
1339 namefull = name.Data();
1340 namefull.Prepend("plots/sigmaVSntracklets_");
1341 namefull.Append(".ps");
1342 c6->Print(namefull.Data());
1343 namefull = name.Data();
1344 namefull.Prepend("plots/sigmaVSntracklets_");
1345 namefull.Append(".eps");
1346 c6->Print(namefull.Data());
1347 namefull = name.Data();
1348 namefull.Prepend("plots/sigmaVSntracklets_");
1349 namefull.Append(".root");
1350 c6->Print(namefull.Data());
1353 namefull = name.Data();
1354 namefull.Prepend("plots/pullsVSntracklets_");
1355 namefull.Append(".gif");
1356 c8->Print(namefull.Data());
1357 namefull = name.Data();
1358 namefull.Prepend("plots/pullsVSntracklets_");
1359 namefull.Append(".ps");
1360 c8->Print(namefull.Data());
1361 namefull = name.Data();
1362 namefull.Prepend("plots/pullsVSntracklets_");
1363 namefull.Append(".eps");
1364 c8->Print(namefull.Data());
1365 namefull = name.Data();
1366 namefull.Prepend("plots/pullsVSntracklets_");
1367 namefull.Append(".root");
1368 c8->Print(namefull.Data());
1373 //----------------------------------------------------------------------------
1374 void TestRmTracks(Int_t iev=0) {
1377 delete gAlice->GetRunLoader();
1382 AliRunLoader *rl = AliRunLoader::Open("galice.root");
1384 cerr<<"Can not open session"<<endl;
1387 Int_t retval = rl->LoadgAlice();
1389 cerr<<"LoadgAlice returned error"<<endl;
1393 retval = rl->LoadHeader();
1395 cerr<<"LoadHeader returned error"<<endl;
1399 gAlice=rl->GetAliRun();
1400 // Get field from galice.root
1401 AliMagF *fiel = (AliMagF*)gAlice->Field();
1402 // Set the conversion constant between curvature and Pt
1403 AliTracker::SetFieldMap(fiel,kTRUE);
1405 Double_t truevtx[3];
1408 Double_t diffX,diffY,diffZ;
1409 Double_t diffXerr,diffYerr,diffZerr;
1410 Double_t multiplicity;
1413 // ----------- Create vertexer --------------------------
1414 AliESDVertex *initVertex = 0;
1415 TFile *invtx = new TFile("AliESDVertexMean.root");
1416 initVertex = (AliESDVertex*)invtx->Get("vtxmean");
1420 AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
1421 vertexer->SetVtxStart(initVertex);
1422 vertexer->SetOnlyFitter();
1423 vertexer->SetDebug(0); // set to 1 to see what it does
1425 Float_t diamondxy[2];
1426 diamondxy[0]=initVertex->GetXv();
1427 diamondxy[1]=initVertex->GetYv();
1428 //----------------------------------------------------------
1430 if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
1431 printf("AliESDs.root not found!\n");
1434 TFile *fin = TFile::Open("AliESDs.root");
1435 AliESDEvent *esd = new AliESDEvent();
1436 TTree *esdTree = (TTree*)fin->Get("esdTree");
1437 if(!esdTree) return;
1438 esd->ReadFromTree(esdTree);
1445 esdTree->GetEvent(iev);
1446 Int_t ntracks = esd->GetNumberOfTracks();
1447 printf("Number of tracks in ESD %d : %d\n",iev,ntracks);
1449 // example: do vertex without tracks 1 and 2 (AliESDtrack::GetID())
1451 // 1: tell vertexer to skip those two tracks
1457 vertexer->SetSkipTracks(nskip,skipped);
1458 AliESDVertex *vertex1 = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
1459 vertex1->PrintStatus();
1460 vertex1->PrintIndices();
1464 // 2: do vertex with all tracks
1468 vertexer->SetSkipTracks(nskip,skipped);
1469 AliESDVertex *vertex2 = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
1470 vertex2->PrintStatus();
1471 vertex2->PrintIndices();
1473 // 3: now remove those two tracks from vertex2
1475 TTree *trksTree = new TTree("trksTree","tree with tracks to be removed");
1476 AliESDtrack *esdTrack = 0;
1477 trksTree->Branch("tracks","AliESDtrack",&esdTrack);
1478 AliESDtrack *et = esd->GetTrack(1);
1479 esdTrack = new AliESDtrack(*et);
1482 AliESDtrack *et = esd->GetTrack(2);
1483 esdTrack = new AliESDtrack(*et);
1486 AliVertexerTracks vt(AliTracker::GetBz());
1487 AliESDVertex *vertex3 = vt.RemoveTracksFromVertex(vertex2,trksTree,diamondxy);
1488 vertex3->PrintStatus();
1489 vertex3->PrintIndices();
1504 //----------------------------------------------------------------------------
1505 void AliComputeVtxMeanFromESD(TString file="AliESDs.root",
1506 Int_t mincontr=20) {
1507 //-----------------------------------------------------------------------
1508 // Compute weighted mean and cov. matrix from a set of AliESD
1509 //-----------------------------------------------------------------------
1511 gStyle->SetOptStat(0);
1512 //gStyle->SetOptFit(0);
1515 Double_t vtx[3],covvtx[6];
1517 AliESDEvent *esd = new AliESDEvent();
1518 AliESDVertex *vertex = 0;
1520 Int_t nc,events,total=0;
1529 Double_t sumwgtx=0.;
1530 Double_t sumwgty=0.;
1531 Double_t sumwgtz=0.;
1541 Double_t eavx,eavy,eavz,ewgtavx,ewgtavy,ewgtavz;
1543 TH1F* hx = new TH1F("hx","",200,-1000,1000);
1544 hx->SetXTitle("vertex x [#mu m]");
1545 hx->SetYTitle("events");
1546 TH1F* hy = new TH1F("hy","",200,-1000,1000);
1547 hy->SetXTitle("vertex y [#mu m]");
1548 hy->SetYTitle("events");
1549 TH1F* hz = new TH1F("hz","",200,-10,10);
1550 hz->SetXTitle("vertex z [cm]");
1551 hz->SetYTitle("events");
1554 for(Int_t i=1; i<100; i++) {
1555 TString inname="event.";
1558 inname.Append(file.Data());
1559 if(gSystem->AccessPathName(inname.Data())) continue;
1560 TFile *fin = TFile::Open(inname.Data());
1561 esdTree = (TTree*)fin->Get("esdTree");
1562 if(!esdTree) continue;
1563 esd->ReadFromTree(esdTree);
1564 events = esdTree->GetEntries();
1565 printf("Number of events in ESD tree: %d\n",events);
1568 for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
1569 esdTree->GetEvent(iev);
1570 vertex = (AliESDVertex*)esd->GetPrimaryVertex();
1571 if(!vertex) continue;
1572 nc = vertex->GetNContributors();
1573 if(nc<mincontr) continue;
1574 vertex->GetXYZ(vtx);
1575 vertex->GetSigmaXYZ(errvtx);
1576 //printf("%f %f %f %f %f %f\n",vtx[0],vtx[1],vtx[2],errvtx[0],errvtx[1],errvtx[2]);
1582 wgtavx += vtx[0]/errvtx[0]/errvtx[0];
1583 wgtavy += vtx[1]/errvtx[1]/errvtx[1];
1584 wgtavz += vtx[2]/errvtx[2]/errvtx[2];
1585 sumwgtx += 1./errvtx[0]/errvtx[0];
1586 sumwgty += 1./errvtx[1]/errvtx[1];
1587 sumwgtz += 1./errvtx[2]/errvtx[2];
1589 hx->Fill(1e4*vtx[0]);
1590 hy->Fill(1e4*vtx[1]);
1606 ewgtavx = 1./TMath::Sqrt(sumwgtx);
1607 ewgtavy = 1./TMath::Sqrt(sumwgty);
1608 ewgtavz = 1./TMath::Sqrt(sumwgtz);
1611 for(Int_t i=1; i<100; i++) {
1612 TString inname="event.";
1615 inname.Append(file.Data());
1616 if(gSystem->AccessPathName(inname.Data())) continue;
1617 TFile *fin = TFile::Open(inname.Data());
1618 esdTree = (TTree*)fin->Get("esdTree");
1619 if(!esdTree) continue;
1620 esd->ReadFromTree(esdTree);
1621 events = esdTree->GetEntries();
1622 for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
1623 esdTree->GetEvent(iev);
1624 vertex = (AliESDVertex*)esd->GetPrimaryVertex();
1625 if(!vertex) continue;
1626 nc = vertex->GetNContributors();
1627 if(nc<mincontr) continue;
1628 vertex->GetXYZ(vtx);
1629 vertex->GetSigmaXYZ(errvtx);
1631 varx += (vtx[0]-avx)*(vtx[0]-avx);
1632 vary += (vtx[1]-avy)*(vtx[1]-avy);
1633 varz += (vtx[2]-avz)*(vtx[2]-avz);
1634 covxy += (vtx[0]-avx)*(vtx[1]-avy);
1635 covxz += (vtx[0]-avx)*(vtx[2]-avz);
1636 covyz += (vtx[1]-avy)*(vtx[2]-avz);
1648 rmsx = TMath::Sqrt(varx);
1649 rmsy = TMath::Sqrt(vary);
1650 rmsz = TMath::Sqrt(varz);
1651 eavx = rmsx/TMath::Sqrt(sum);
1652 eavy = rmsy/TMath::Sqrt(sum);
1653 eavz = rmsz/TMath::Sqrt(sum);
1656 printf("\n\nNumber of events: Total %d, Used %d\n",total,sum);
1657 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);
1658 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);
1659 printf("RMS:\n x = %f mum\n y = %f mum\n z = %f mum\n",1e4*rmsx,1e4*rmsy,1e4*rmsz);
1661 TCanvas *c = new TCanvas("c","c",0,0,1000,500);
1665 TF1 *gx = new TF1("gx","gaus",-1000,1000);
1666 gx->SetLineColor(2);
1668 TF1 *gxx = (TF1*)gx->Clone("gxx");
1669 gxx->FixParameter(2,50.);
1670 gxx->SetLineStyle(2);
1674 TF1 *gy = new TF1("gy","gaus",-1000,1000);
1675 gy->SetLineColor(2);
1677 TF1 *gyy = (TF1*)gy->Clone("gyy");
1678 gyy->FixParameter(2,50.);
1679 gyy->SetLineStyle(2);
1683 TF1 *gz = new TF1("gz","gaus",-10,10);
1684 gz->SetLineColor(2);
1686 TF1 *gzz = (TF1*)gz->Clone("gzz");
1687 gzz->FixParameter(2,5.3);
1688 gzz->SetLineStyle(2);
1700 AliESDVertex *vtxmean = new AliESDVertex(vtx,covvtx,0.,(Int_t)sum);
1701 vtxmean->SetTitle("Vertex from weighted average");
1702 TFile *out = new TFile("AliESDVertexMean.root","recreate");
1703 vtxmean->Write("vtxmean");
1708 //----------------------------------------------------------------------------
1709 void AliComputeVtxMeanFromTree(TString file="AliESDVertices.root",
1710 Int_t mincontr=20) {
1711 //-----------------------------------------------------------------------
1712 // Compute weighted mean and cov. matrix from a set of AliESDVertices
1713 //-----------------------------------------------------------------------
1714 gStyle->SetOptStat(0);
1715 //gStyle->SetOptFit(0);
1718 Double_t vtx[3],covvtx[6];
1729 Double_t sumwgtx=0.;
1730 Double_t sumwgty=0.;
1731 Double_t sumwgtz=0.;
1735 Double_t eavx,eavy,eavz,ewgtavx,ewgtavy,ewgtavz;
1737 TH1F* hx = new TH1F("hx","",200,-1000,1000);
1738 hx->SetXTitle("vertex x [#mu m]");
1739 hx->SetYTitle("events");
1740 TH1F* hy = new TH1F("hy","",200,-1000,1000);
1741 hy->SetXTitle("vertex y [#mu m]");
1742 hy->SetYTitle("events");
1743 TH1F* hz = new TH1F("hz","",200,-10,10);
1744 hz->SetXTitle("vertex z [cm]");
1745 hz->SetYTitle("events");
1748 if(gSystem->AccessPathName(file.Data())) continue;
1749 TFile *fin = TFile::Open(file.Data());
1750 TTree *tree = (TTree*)fin->Get("TreeV");
1751 AliESDVertex *vertex = 0;
1752 tree->SetBranchAddress("AliESDVertex",&vertex);
1753 Int_t entries = (Int_t)tree->GetEntries();
1755 for(Int_t i=0; i<entries; i++) {
1758 nc = vertex->GetNContributors();
1759 if(nc<mincontr) continue;
1760 vertex->GetXYZ(vtx);
1761 vertex->GetSigmaXYZ(errvtx);
1762 //printf("%f %f %f %f %f %f\n",vtx[0],vtx[1],vtx[2],errvtx[0],errvtx[1],errvtx[2]);
1768 wgtavx += vtx[0]/errvtx[0]/errvtx[0];
1769 wgtavy += vtx[1]/errvtx[1]/errvtx[1];
1770 wgtavz += vtx[2]/errvtx[2]/errvtx[2];
1771 sumwgtx += 1./errvtx[0]/errvtx[0];
1772 sumwgty += 1./errvtx[1]/errvtx[1];
1773 sumwgtz += 1./errvtx[2]/errvtx[2];
1775 hx->Fill(1e4*vtx[0]);
1776 hy->Fill(1e4*vtx[1]);
1788 ewgtavx = 1./TMath::Sqrt(sumwgtx);
1789 ewgtavy = 1./TMath::Sqrt(sumwgty);
1790 ewgtavz = 1./TMath::Sqrt(sumwgtz);
1793 for(Int_t i=0; i<entries; i++) {
1795 nc = vertex->GetNContributors();
1796 if(nc<mincontr) continue;
1797 vertex->GetXYZ(vtx);
1798 vertex->GetSigmaXYZ(errvtx);
1800 varx += (vtx[0]-avx)*(vtx[0]-avx);
1801 vary += (vtx[1]-avy)*(vtx[1]-avy);
1802 varz += (vtx[2]-avz)*(vtx[2]-avz);
1803 covxy += (vtx[0]-avx)*(vtx[1]-avy);
1804 covxz += (vtx[0]-avx)*(vtx[2]-avz);
1805 covyz += (vtx[1]-avy)*(vtx[2]-avz);
1816 rmsx = TMath::Sqrt(varx);
1817 rmsy = TMath::Sqrt(vary);
1818 rmsz = TMath::Sqrt(varz);
1819 eavx = rmsx/TMath::Sqrt(sum);
1820 eavy = rmsy/TMath::Sqrt(sum);
1821 eavz = rmsz/TMath::Sqrt(sum);
1824 printf("\n\nNumber of events: Total %d, Used %d\n",total,sum);
1825 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);
1826 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);
1827 printf("RMS:\n x = %f mum\n y = %f mum\n z = %f mum\n",1e4*rmsx,1e4*rmsy,1e4*rmsz);
1829 TCanvas *c = new TCanvas("c","c",0,0,1000,500);
1833 TF1 *gx = new TF1("gx","gaus",-1000,1000);
1834 gx->SetLineColor(2);
1836 TF1 *gxx = (TF1*)gx->Clone("gxx");
1837 gxx->FixParameter(2,50.);
1838 gxx->SetLineStyle(2);
1842 TF1 *gy = new TF1("gy","gaus",-1000,1000);
1843 gy->SetLineColor(2);
1845 TF1 *gyy = (TF1*)gy->Clone("gyy");
1846 gyy->FixParameter(2,50.);
1847 gyy->SetLineStyle(2);
1851 TF1 *gz = new TF1("gz","gaus",-10,10);
1852 gz->SetLineColor(2);
1854 TF1 *gzz = (TF1*)gz->Clone("gzz");
1855 gzz->FixParameter(2,5.3);
1856 gzz->SetLineStyle(2);
1868 AliESDVertex *vtxmean = new AliESDVertex(vtx,covvtx,0.,(Int_t)sum);
1869 vtxmean->SetTitle("Vertex from weighted average");
1870 TFile *out = new TFile("AliESDVertexMean.root","recreate");
1871 vtxmean->Write("vtxmean");
1876 //---------------------------------------------------------------------------