ALIROOT-5486 Propagate MUON trigger deviation sign to AliAODTrack
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracksTest.C
CommitLineData
43c9dae1 1void AliVertexerTracksTest(TString outname="AliVertexerTracksTest.root",
2 Bool_t tpconly=kTRUE,
fee9ef0d 3 Bool_t useMeanVtx=kFALSE,
43c9dae1 4 Double_t nSigma=3.,
fee9ef0d 5 Int_t itsrefit=1,
6 Double_t maxd0z0=0.5,
7 Int_t minitscls=5,
8 Int_t mintracks=1,
fee9ef0d 9 Bool_t onlyfit=kFALSE) {
820252d7 10//------------------------------------------------------------------------
11// Run vertex reconstruction and store results in histos and ntuple
12//------------------------------------------------------------------------
fee9ef0d 13// WITHOUT Kinematics.root
14//
43c9dae1 15 AliGeomManager::LoadGeometry("geometry.root");
16
820252d7 17 if (gAlice) {
33c3c91a 18 delete AliRunLoader::Instance();
820252d7 19 delete gAlice;
20 gAlice=0;
21 }
22
23 AliRunLoader *rl = AliRunLoader::Open("galice.root");
24 if (rl == 0x0) {
25 cerr<<"Can not open session"<<endl;
26 return;
27 }
28 Int_t retval = rl->LoadgAlice();
29 if (retval) {
30 cerr<<"LoadgAlice returned error"<<endl;
31 delete rl;
32 return;
33 }
34 retval = rl->LoadHeader();
35 if (retval) {
36 cerr<<"LoadHeader returned error"<<endl;
37 delete rl;
38 return;
39 }
40 gAlice=rl->GetAliRun();
fee9ef0d 41
42 //rl->LoadKinematics();
43
820252d7 44 // Get field from galice.root
45 AliMagF *fiel = (AliMagF*)gAlice->Field();
46 // Set the conversion constant between curvature and Pt
47 AliTracker::SetFieldMap(fiel,kTRUE);
48
49 //**** Switch on the PID class
50 AliPID pid;
51
52 Int_t nEvents = (Int_t)rl->GetNumberOfEvents();
53
54 Double_t truevtx[3];
55 Double_t vtx[3];
56 Double_t errvtx[3];
57 Double_t diffX,diffY,diffZ;
58 Double_t diffXerr,diffYerr,diffZerr;
fee9ef0d 59 Float_t multiplicity;
60 Float_t chi2red;
61 Int_t ntracklets;
820252d7 62 Int_t nc;
fee9ef0d 63 Int_t triggered;
820252d7 64
65 // Histograms
66 TH1F *hdiffX = new TH1F("hdiffX","x_{found} - x_{true} [#mum]",1000,-5000,5000);
67 TH1F *hdiffY = new TH1F("hdiffY","y_{found} - y_{true} [#mum]",1000,-5000,5000);
68 TH1F *hdiffZ = new TH1F("hdiffZ","z_{found} - z_{true} [#mum]",1000,-5000,5000);
69
70 TH1F *hdiffXerr = new TH1F("hdiffXerr","#Delta x/#sigma_{x}",100,-5,5);
71 TH1F *hdiffYerr = new TH1F("hdiffYerr","#Delta y/#sigma_{y}",100,-5,5);
72 TH1F *hdiffZerr = new TH1F("hdiffZerr","#Delta z/#sigma_{z}",100,-5,5);
73
74 //TNtple
fee9ef0d 75 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");
820252d7 76
77
78
79 // ----------- Create vertexer --------------------------
80 AliESDVertex *initVertex = 0;
81 if(useMeanVtx) {
82 // open the mean vertex
83 TFile *invtx = new TFile("AliESDVertexMean.root");
84 initVertex = (AliESDVertex*)invtx->Get("vtxmean");
85 invtx->Close();
86 delete invtx;
87 } else {
88 Double_t pos[3]={0.5,0.5,0.};
89 Double_t err[3]={3.,3.,5.};
90 initVertex = new AliESDVertex(pos,err);
91 }
28e4388a 92 AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
fee9ef0d 93 vertexer->SetDebug(1); // set to 1 to see what it does
820252d7 94 vertexer->SetVtxStart(initVertex);
0bf2e2b4 95 if(!useMeanVtx) vertexer->SetConstraintOff();
fee9ef0d 96 if(onlyfit) vertexer->SetOnlyFitter();
820252d7 97 vertexer->SetNSigmad0(nSigma);
43c9dae1 98 if(!itsrefit || tpconly) vertexer->SetITSrefitNotRequired();
fee9ef0d 99 vertexer->SetMinTracks(mintracks);
100 vertexer->SetMinITSClusters(minitscls);
101 vertexer->SetMaxd0z0(maxd0z0);
820252d7 102 //cout<<" Nsigma: "<<vertexer->GetNSigmad0()<<endl;
820252d7 103 //----------------------------------------------------------
104
fee9ef0d 105
43c9dae1 106 if(gSystem->AccessPathName("AliESDs_itstpc.root",kFileExists)) {
107 printf("AliESDs_itstpc.root not found!\n");
820252d7 108 return;
109 }
43c9dae1 110 TFile *fin = TFile::Open("AliESDs_itstpc.root");
64cf164e 111 AliESDEvent *esd = new AliESDEvent();
820252d7 112 TTree *esdTree = (TTree*)fin->Get("esdTree");
113 if(!esdTree) return;
64cf164e 114 esd->ReadFromTree(esdTree);
820252d7 115 Int_t events = esdTree->GetEntries();
116 printf("Number of events in ESD tree: %d\n",events);
117
118 TArrayF o;
119
120 for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
121 printf("-------------- EVENT %d --------------------\n",iev);
122
43c9dae1 123 diffX=999999;
124 diffY=999999;
125 diffZ=999999;
126 diffXerr=999999;
127 diffYerr=999999;
128 diffZerr=999999;
820252d7 129 nc=0;
130
131 //=================================================
132 // LOOK IN THE SIMULATION
133 // true vertex position
134 Int_t npart = (Int_t)gAlice->GetEvent(iev);
135 printf("particles %d\n",npart);
136 AliHeader *header = (AliHeader*)rl->GetHeader();
137 AliGenPythiaEventHeader *pyheader = (AliGenPythiaEventHeader*)header->GenEventHeader();
138 pyheader->PrimaryVertex(o);
139 truevtx[0] = o[0];
140 truevtx[1] = o[1];
141 truevtx[2] = o[2];
142 printf("True Vertex: (%f, %f, %f)\n",o[0],o[1],o[2]);
143
144 Double_t dNchdy = 0.;
145 multiplicity=0.;
146 Int_t nitstracks = 0;
820252d7 147 Int_t nitstracks1 = 0;
148 Int_t nitstracks2 = 0;
149 Int_t nitstracks3 = 0;
150 Int_t nitstracks4 = 0;
151 Int_t nitstracks5 = 0;
152 Int_t nitstracks6 = 0;
153 Int_t nitstracks5or6 = 0;
820252d7 154
155
156 esdTree->GetEvent(iev);
fee9ef0d 157 triggered=0;
158 chi2red=0.;
159 ULong64_t triggerMask = esd->GetTriggerMask();
160 if(triggerMask&32 && ((triggerMask&1) || (triggerMask&2))) triggered=1;
161
162 // get # of tracklets
163 AliMultiplicity *alimult = esd->GetMultiplicity();
164 if(alimult) {
165 ntracklets = alimult->GetNumberOfTracklets();
166 for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++)
167 if(alimult->GetDeltaPhi(l)<-9998.) ntracklets--;
168 } else {
169 ntracklets = 0;
170 }
171 printf("Number of SPD tracklets in ESD %d : %d\n",iev,ntracklets);
172 multiplicity = (Float_t)ntracklets;
173
820252d7 174 Int_t ntracks = esd->GetNumberOfTracks();
175 printf("Number of tracks in ESD %d : %d\n",iev,ntracks);
176 if(ntracks==0) {
fee9ef0d 177 ntVtxRes->Fill(triggered,ntracks,nitstracks,nitstracks5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,nc,truevtx[2],chi2red);
820252d7 178 continue;
179 }
180
181 // VERTEX
43c9dae1 182 AliESDVertex *vertex = 0;
183 if(!tpconly) {
184 vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
185 } else {
186 TObjArray trkArrayOrig(ntracks);
187 UShort_t *idOrig = new UShort_t[ntracks];
188 const Double_t kRadius = 2.8; //something less than the beam pipe radius
189 const Double_t kMaxStep = 5; //max step over the material
190 Bool_t ok;
191 Int_t nTrksOrig=0;
192 for(Int_t i=0; i<ntracks; i++) {
193 AliESDtrack *esdt = esd->GetTrack(i);
194 AliExternalTrackParam *tpcTrack =
195 (AliExternalTrackParam *)esdt->GetTPCInnerParam();
196 ok = kFALSE;
197 if (tpcTrack)
198 ok = AliTracker::
199 PropagateTrackTo(tpcTrack,kRadius,esdt->GetMass(),kRadius,kTRUE);
200 if (ok) {
201 trkArrayOrig.AddLast(tpcTrack);
202 idOrig[nTrksOrig]=(UShort_t)esdt->GetID();
203 nTrksOrig++;
204 }
205 }
206
207 vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(&trkArrayOrig,idOrig);
208 delete [] idOrig;
209 }
210
820252d7 211
212 nc = (Int_t)vertex->GetNContributors();
fee9ef0d 213 if(nc>=2) chi2red = vertex->GetChi2toNDF();
820252d7 214 printf("nc = %d\n",nc);
215
216
217 // count ITS tracks
218 for(Int_t itrk=0; itrk<ntracks; itrk++) {
219 AliESDtrack *esdtrack = (AliESDtrack*)esd->GetTrack(itrk);
220 UInt_t status = esdtrack->GetStatus();
221 // only tracks found also in ITS
222 if(! (status&AliESDtrack::kITSin) ) continue;
223 nitstracks++;
224 Int_t npts = (Int_t)esdtrack->GetNcls(0);
225 if(npts==6) {nitstracks6++;nitstracks5or6++;}
226 if(npts==5) {nitstracks5++;nitstracks5or6++;}
227 if(npts==4) nitstracks4++;
228 if(npts==3) nitstracks3++;
229 if(npts==2) nitstracks2++;
230 if(npts==1) nitstracks1++;
820252d7 231 }
232 printf("Number of kITSin tracks in ESD %d : %d\n",iev,nitstracks);
233 printf(" 6 pts: %d\n",nitstracks6);
234 printf(" 5 pts: %d\n",nitstracks5);
235 printf(" 4 pts: %d\n",nitstracks4);
236 printf(" 3 pts: %d\n",nitstracks3);
237 printf(" 2 pts: %d\n",nitstracks2);
238 printf(" 1 pts: %d\n",nitstracks1);
820252d7 239
240
fee9ef0d 241 if(nc>=1) {
820252d7 242 vertex->GetXYZ(vtx);
243 vertex->GetSigmaXYZ(errvtx);
244 diffX = 10000.* (vtx[0] - truevtx[0]);
245 diffY = 10000.* (vtx[1] - truevtx[1]);
246 diffZ = 10000.* (vtx[2] - truevtx[2]);
247 hdiffX->Fill(diffX);
248 hdiffY->Fill(diffY);
249 hdiffZ->Fill(diffZ);
250 diffXerr = diffX/errvtx[0]/10000.;
251 diffYerr = diffY/errvtx[1]/10000.;
252 diffZerr = diffZ/errvtx[2]/10000.;
253 hdiffXerr->Fill(diffXerr);
254 hdiffYerr->Fill(diffYerr);
255 hdiffZerr->Fill(diffZerr);
256 }
257
fee9ef0d 258 ntVtxRes->Fill(triggered,ntracks,nitstracks,nitstracks5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,(Float_t)nc,truevtx[2],chi2red);
820252d7 259
260 } // END LOOP ON EVENTS
261
262
263 delete esdTree;
264 delete fin;
265 delete vertexer;
266 delete initVertex;
267
268 // Output file
269 TFile *fout = new TFile(outname.Data(),"recreate");
270 ntVtxRes->Write();
271 hdiffX->Write();
272 hdiffY->Write();
273 hdiffZ->Write();
274 hdiffXerr->Write();
275 hdiffYerr->Write();
276 hdiffZerr->Write();
277 fout->Close();
278
279 return;
280}
fee9ef0d 281//------------------------------------------------------------------------
282void VertexForOneEvent(Int_t iev=0,
283 Double_t nSigma=3.,
284 Bool_t useMeanVtx=kFALSE) {
285//------------------------------------------------------------------------
286// Run vertex reconstruction for event iev
287//------------------------------------------------------------------------
288
289 if (gAlice) {
33c3c91a 290 delete AliRunLoader::Instance();
fee9ef0d 291 delete gAlice;
292 gAlice=0;
293 }
294
295 AliRunLoader *rl = AliRunLoader::Open("galice.root");
296 if (rl == 0x0) {
297 cerr<<"Can not open session"<<endl;
298 return;
299 }
300 Int_t retval = rl->LoadgAlice();
301 if (retval) {
302 cerr<<"LoadgAlice returned error"<<endl;
303 delete rl;
304 return;
305 }
306 retval = rl->LoadHeader();
307 if (retval) {
308 cerr<<"LoadHeader returned error"<<endl;
309 delete rl;
310 return;
311 }
312 gAlice=rl->GetAliRun();
313 //rl->LoadKinematics();
314 // Get field from galice.root
315 AliMagF *fiel = (AliMagF*)gAlice->Field();
316 // Set the conversion constant between curvature and Pt
317 AliTracker::SetFieldMap(fiel,kTRUE);
318
319 Double_t truevtx[3];
320 Double_t vtx[3];
321 Double_t errvtx[3];
322 Double_t diffX,diffY,diffZ;
323 Double_t diffXerr,diffYerr,diffZerr;
324 Double_t multiplicity;
325 Int_t nc;
326
327 // ----------- Create vertexer --------------------------
328 AliESDVertex *initVertex = 0;
329 if(useMeanVtx) {
330 // open the mean vertex
331 TFile *invtx = new TFile("AliESDVertexMean.root");
332 initVertex = (AliESDVertex*)invtx->Get("vtxmean");
333 invtx->Close();
334 delete invtx;
335 } else {
336 Double_t pos[3]={0.,-0.,0.};
337 Double_t err[3]={3.,3.,5.};
338 initVertex = new AliESDVertex(pos,err);
339 }
28e4388a 340 AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
fee9ef0d 341 //vertexer->SetITSrefitNotRequired();
342 vertexer->SetVtxStart(initVertex);
0bf2e2b4 343 if(!useMeanVtx) vertexer->SetConstraintOff();
fee9ef0d 344 //vertexer->SetMinTracks(2);
345 //vertexer->SetNSigmad0(nSigma);
346 //cout<<" Nsigma: "<<vertexer->GetNSigmad0()<<endl;
347 //vertexer->SetFinderAlgorithm(2);
348 //vertexer->SetDCAcut(0.1); // 1 mm
349 vertexer->SetDebug(1); // set to 1 to see what it does
350 //----------------------------------------------------------
351
352 if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
353 printf("AliESDs.root not found!\n");
354 return;
355 }
356 TFile *fin = TFile::Open("AliESDs.root");
64cf164e 357 AliESDEvent *esd = new AliESDEvent();
fee9ef0d 358 TTree *esdTree = (TTree*)fin->Get("esdTree");
359 if(!esdTree) return;
64cf164e 360 esd->ReadFromTree(esdTree);
fee9ef0d 361
362 TArrayF o;
363
364 nc=0;
365
366 //=================================================
367 // LOOK IN THE SIMULATION
368 // true vertex position
369 Int_t npart = (Int_t)gAlice->GetEvent(iev);
370 printf("particles %d\n",npart);
371 AliHeader *header = (AliHeader*)rl->GetHeader();
372 AliGenPythiaEventHeader *pyheader = (AliGenPythiaEventHeader*)header->GenEventHeader();
373 pyheader->PrimaryVertex(o);
374 truevtx[0] = o[0];
375 truevtx[1] = o[1];
376 truevtx[2] = o[2];
377 printf("True Vertex: (%f, %f, %f)\n",o[0],o[1],o[2]);
378
379 esdTree->GetEvent(iev);
380 Int_t ntracks = esd->GetNumberOfTracks();
381 printf("Number of tracks in ESD %d : %d\n",iev,ntracks);
382
383 // VERTEX
384 AliESDVertex *vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
385 nc = (Int_t)vertex->GetNContributors();
386
387
388 Int_t nitstracks = 0;
389 Int_t nitstracks1 = 0;
390 Int_t nitstracks2 = 0;
391 Int_t nitstracks3 = 0;
392 Int_t nitstracks4 = 0;
393 Int_t nitstracks5 = 0;
394 Int_t nitstracks6 = 0;
395 Int_t nitstracks5or6 = 0;
396
397 // count ITS tracks
398 for(Int_t itrk=0; itrk<ntracks; itrk++) {
399 AliESDtrack *esdtrack = (AliESDtrack*)esd->GetTrack(itrk);
400 UInt_t status = esdtrack->GetStatus();
401 // only tracks found also in ITS
402 if(! (status&AliESDtrack::kITSrefit) ) continue;
403 nitstracks++;
404 Int_t npts = (Int_t)esdtrack->GetNcls(0);
405 if(npts==6) {nitstracks6++;nitstracks5or6++;}
406 if(npts==5) {nitstracks5++;nitstracks5or6++;}
407 if(npts==4) nitstracks4++;
408 if(npts==3) nitstracks3++;
409 if(npts==2) nitstracks2++;
410 if(npts==1) nitstracks1++;
411 }
412 printf("Number of kITSrefit tracks in ESD %d : %d\n",iev,nitstracks);
413 printf(" 6 pts: %d\n",nitstracks6);
414 printf(" 5 pts: %d\n",nitstracks5);
415 printf(" 4 pts: %d\n",nitstracks4);
416 printf(" 3 pts: %d\n",nitstracks3);
417 printf(" 2 pts: %d\n",nitstracks2);
418 printf(" 1 pts: %d\n",nitstracks1);
419
420
421 delete esdTree;
422 delete fin;
423 delete vertexer;
424 delete initVertex;
425
426 return;
427}
428//----------------------------------------------------------------------------
820252d7 429Double_t FitFunc(Double_t *x,Double_t *par) {
430 return par[0]+par[1]/TMath::Sqrt(x[0]);
431}
432Int_t GetBin(Float_t mult) {
433 if(mult<5.) return 0;
434 if(mult<7.) return 1;
435 if(mult<12.) return 2;
436 if(mult<15.) return 3;
437 if(mult<22.) return 4;
438 return 5;
439}
fee9ef0d 440Int_t GetBinTracklets(Float_t tracklets) {
441 if(tracklets<2.*5.) return 0;
442 if(tracklets<2.*7.) return 1;
443 if(tracklets<2.*12.) return 2;
444 if(tracklets<2.*15.) return 3;
445 if(tracklets<2.*22.) return 4;
446 return 5;
447}
448Int_t GetBinZ(Float_t z) {
449 if(z<-13.) return 0;
450 if(z<-11.) return 1;
451 if(z<-9.) return 2;
452 if(z<-7.) return 3;
453 if(z<-5.) return 4;
454 if(z<-3.) return 5;
455 if(z<-1.) return 6;
456 if(z<1.) return 7;
457 if(z<3.) return 8;
458 if(z<5.) return 9;
459 if(z<7.) return 10;
460 if(z<9.) return 11;
461 if(z<11.) return 12;
462 if(z<13.) return 13;
463 return 14;
464}
465//----------------------------------------------------------------------------
43c9dae1 466void PlotVtxRes(TString inname="AliVertexerTracksTest.root",
467 Bool_t tpconly=kFALSE) {
820252d7 468 //---------------------------------------------------------------------
469 // Plot efficiency, resolutions, pulls
470 // [at least 10k pp events are needed]
471 //---------------------------------------------------------------------
472
473 gStyle->SetOptStat(0);
474 gStyle->SetOptFit(10001);
475
43c9dae1 476 Float_t zMCmin=0.;
477 Float_t zMCmax=15.;
478 Float_t ncminforzshift=1.;
479 Float_t maxdiffx = 1e3;
480 Float_t rangefactor=1.;
481 if(tpconly) {
482 maxdiffx *= 1e2;
483 rangefactor = 15.;
484 }
485
820252d7 486 Int_t nEvVtx=0;
487 Int_t nEvNoVtx=0;
488 Int_t ev,nUsedTrks;
fee9ef0d 489 Float_t nESDTrks,nTrks,nTrks5or6,ntracklets;
490 Float_t diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,zMC,nc;
491 Float_t triggered;
492
820252d7 493
494 //
495 // HISTOGRAMS
496 //
43c9dae1 497 TH1F *hdx = new TH1F("hdx","",50,-600*rangefactor,600*rangefactor);
820252d7 498 hdx->SetXTitle("#Delta x [#mu m]");
499 hdx->SetYTitle("events");
500 //
43c9dae1 501 TH1F *hdy = new TH1F("hdy","",50,-600*rangefactor,600*rangefactor);
820252d7 502 hdy->SetXTitle("#Delta y [#mu m]");
503 hdy->SetYTitle("events");
504 //
43c9dae1 505 TH1F *hdz = new TH1F("hdz","",200,-600*rangefactor,600*rangefactor);
820252d7 506 hdz->SetXTitle("#Delta z [#mu m]");
507 hdz->SetYTitle("events");
508 //
fee9ef0d 509 TH1F *hzMCVtx = new TH1F("hzMCVtx","events w/ vertex",30,-15,15);
510 hzMCVtx->SetXTitle("z_{true} [cm]");
511 hzMCVtx->SetYTitle("events");
512 hzMCVtx->Sumw2();
513 //
514 TH1F *hzMCNoVtx = new TH1F("hzMCNoVtx","events w/o vertex",30,-15,15);
515 hzMCNoVtx->SetXTitle("z_{true} [cm]");
516 hzMCNoVtx->SetYTitle("events");
517 hzMCNoVtx->Sumw2();
518 //
519 TH1F *hTrackletsVtx = new TH1F("hTrackletsVtx","events w/ vertex",51,-0.5,50.5);
520 hTrackletsVtx->SetXTitle("SPD tracklets");
521 hTrackletsVtx->SetYTitle("events");
522 //
523 TH1F *hTrackletsNoVtx = new TH1F("hTrackletsNoVtx","events w/o vertex",51,-0.5,50.5);
524 hTrackletsNoVtx->SetXTitle("SPD tracklets");
525 hTrackletsNoVtx->SetYTitle("events");
526 //
820252d7 527 TH1F *hTracksVtx = new TH1F("hTracksVtx","events w/ vertex",51,-0.5,50.5);
528 hTracksVtx->SetXTitle("Number of reco tracks in TPC&ITS");
529 hTracksVtx->SetYTitle("events");
530 //
531 TH1F *hTracksNoVtx = new TH1F("hTracksNoVtx","events w/o vertex",51,-0.5,50.5);
532 hTracksNoVtx->SetXTitle("Number of reco tracks in TPC&ITS");
533 hTracksNoVtx->SetYTitle("events");
534 //
535 TH1F *hTracksVtx5or6 = new TH1F("hTracksVtx5or6","events w/ vertex",51,-0.5,50.5);
fee9ef0d 536 hTracksVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(cls>=5)");
820252d7 537 hTracksVtx5or6->SetYTitle("events");
538 //
539 TH1F *hTracksNoVtx5or6 = new TH1F("hTracksNoVtx5or6","events w/o vertex",51,-0.5,50.5);
fee9ef0d 540 hTracksNoVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(cls>=5)");
820252d7 541 hTracksNoVtx5or6->SetYTitle("events");
542 //
820252d7 543 TH1F *hPullsx = new TH1F("hPullsx","",50,-7.,7.);
544 hPullsx->SetXTitle("#Delta x/#sqrt{<x,x>}");
545 hPullsx->SetYTitle("events");
546 //
547 TH1F *hPullsy = new TH1F("hPullsy","",50,-7.,7.);
548 hPullsy->SetXTitle("#Delta y/#sqrt{<y,y>}");
549 hPullsy->SetYTitle("events");
550 //
551 TH1F *hPullsz = new TH1F("hPullsz","",50,-7.,7.);
552 hPullsz->SetXTitle("#Delta z/#sqrt{<z,z>}");
553 hPullsz->SetYTitle("events");
554 //
fee9ef0d 555 TProfile *hntrackletsVSzMC = new TProfile("hntrackletsVSzMC","",30,-15,15,0,30);
556 hntrackletsVSzMC->SetXTitle("z_{true} [cm]");
557 hntrackletsVSzMC->SetYTitle("<SPD tracklets>");
558 //
559 TProfile *hnitstracksVSzMC = new TProfile("hnitstracksVSzMC","",30,-15,15,0,30);
560 hnitstracksVSzMC->SetXTitle("z_{true} [cm]");
561 hnitstracksVSzMC->SetYTitle("<tracks TPC&ITS>");
562 //
563 TProfile *hnitstracks5or6VSzMC = new TProfile("hnitstracks5or6VSzMC","",30,-15,15,0,30);
564 hnitstracks5or6VSzMC->SetXTitle("z_{true} [cm]");
565 hnitstracks5or6VSzMC->SetYTitle("<tracks TPC&ITS(cls>=5)>");
820252d7 566
567 Double_t mult[6]={0,0,0,0,0,0};
568 Double_t mult2[6]={0,0,0,0,0,0};
569 Double_t sigmamult[6]={0,0,0,0,0,0};
570 Double_t sigmamult2[6]={0,0,0,0,0,0};
571 Double_t tottracks[6]={0,0,0,0,0,0};
572 Double_t fittrks[6]={0,0,0,0,0,0};
573 Double_t tracks[6]={0,0,0,0,0,0};
574 Double_t etracks[6]={0,0,0,0,0,0};
575 Double_t xres[6]={0,0,0,0,0,0};
576 Double_t exres[6]={0,0,0,0,0,0};
577 Double_t yres[6]={0,0,0,0,0,0};
578 Double_t eyres[6]={0,0,0,0,0,0};
579 Double_t zres[6]={0,0,0,0,0,0};
580 Double_t ezres[6]={0,0,0,0,0,0};
581 Double_t xpull[6]={0,0,0,0,0,0};
582 Double_t expull[6]={0,0,0,0,0,0};
583 Double_t ypull[6]={0,0,0,0,0,0};
584 Double_t eypull[6]={0,0,0,0,0,0};
585 Double_t zpull[6]={0,0,0,0,0,0};
586 Double_t ezpull[6]={0,0,0,0,0,0};
587 Int_t evts[6]={0,0,0,0,0,0};
588 Int_t totevts[6]={0,0,0,0,0,0};
589 Int_t vtx[6]={0,0,0,0,0,0};
590 Int_t all[6]={0,0,0,0,0,0};
591 Double_t probVtx[6]={0,0,0,0,0,0};
592 Double_t eprobVtx[6]={0,0,0,0,0,0};
fee9ef0d 593 Int_t bin,zbin;
820252d7 594 //
595 // OTHER HISTOGRAMS
596 //
43c9dae1 597 TH1F *hdx0 = new TH1F("hdx0","x resolution - bin 0",50,-500*rangefactor,500*rangefactor);
598 TH1F *hdx1 = new TH1F("hdx1","x resolution - bin 1",50,-500*rangefactor,500*rangefactor);
599 TH1F *hdx2 = new TH1F("hdx2","x resolution - bin 2",50,-500*rangefactor,500*rangefactor);
600 TH1F *hdx3 = new TH1F("hdx3","x resolution - bin 3",50,-400*rangefactor,400*rangefactor);
601 TH1F *hdx4 = new TH1F("hdx4","x resolution - bin 4",50,-300*rangefactor,300*rangefactor);
602 TH1F *hdx5 = new TH1F("hdx5","x resolution - bin 5",50,-300*rangefactor,300*rangefactor);
603 TH1F *hdy0 = new TH1F("hdy0","y resolution - bin 0",50,-500*rangefactor,500*rangefactor);
604 TH1F *hdy1 = new TH1F("hdy1","y resolution - bin 1",50,-500*rangefactor,500*rangefactor);
605 TH1F *hdy2 = new TH1F("hdy2","y resolution - bin 2",50,-500*rangefactor,500*rangefactor);
606 TH1F *hdy3 = new TH1F("hdy3","y resolution - bin 3",50,-400*rangefactor,400*rangefactor);
607 TH1F *hdy4 = new TH1F("hdy4","y resolution - bin 4",50,-300*rangefactor,300*rangefactor);
608 TH1F *hdy5 = new TH1F("hdy5","y resolution - bin 5",50,-300*rangefactor,300*rangefactor);
609 TH1F *hdz0 = new TH1F("hdz0","z resolution - bin 0",50,-1000*rangefactor,1000*rangefactor);
610 TH1F *hdz1 = new TH1F("hdz1","z resolution - bin 1",50,-800*rangefactor,800*rangefactor);
611 TH1F *hdz2 = new TH1F("hdz2","z resolution - bin 2",50,-800*rangefactor,800*rangefactor);
612 TH1F *hdz3 = new TH1F("hdz3","z resolution - bin 3",50,-600*rangefactor,600*rangefactor);
613 TH1F *hdz4 = new TH1F("hdz4","z resolution - bin 4",50,-500*rangefactor,500*rangefactor);
614 TH1F *hdz5 = new TH1F("hdz5","z resolution - bin 5",50,-500*rangefactor,500*rangefactor);
fee9ef0d 615 //
616 TH1F *hdz_z = NULL; hdz_z = new TH1F[15];
617 for(Int_t l=0;l<15;l++) hdz_z[l] = *hdz;
618
820252d7 619
620 TH1F *hpx0 = new TH1F("hpx0","x pull - bin 0",50,-7,7);
621 TH1F *hpx1 = new TH1F("hpx1","x pull - bin 1",50,-7,7);
622 TH1F *hpx2 = new TH1F("hpx2","x pull - bin 2",50,-7,7);
623 TH1F *hpx3 = new TH1F("hpx3","x pull - bin 3",50,-7,7);
624 TH1F *hpx4 = new TH1F("hpx4","x pull - bin 4",50,-7,7);
625 TH1F *hpx5 = new TH1F("hpx5","x pull - bin 5",50,-7,7);
626 TH1F *hpy0 = new TH1F("hpy0","y pull - bin 0",50,-7,7);
627 TH1F *hpy1 = new TH1F("hpy1","y pull - bin 1",50,-7,7);
628 TH1F *hpy2 = new TH1F("hpy2","y pull - bin 2",50,-7,7);
629 TH1F *hpy3 = new TH1F("hpy3","y pull - bin 3",50,-7,7);
630 TH1F *hpy4 = new TH1F("hpy4","y pull - bin 4",50,-7,7);
631 TH1F *hpy5 = new TH1F("hpy5","y pull - bin 5",50,-7,7);
632 TH1F *hpz0 = new TH1F("hpz0","z pull - bin 0",50,-7,7);
633 TH1F *hpz1 = new TH1F("hpz1","z pull - bin 1",50,-7,7);
634 TH1F *hpz2 = new TH1F("hpz2","z pull - bin 2",50,-7,7);
635 TH1F *hpz3 = new TH1F("hpz3","z pull - bin 3",50,-7,7);
636 TH1F *hpz4 = new TH1F("hpz4","z pull - bin 4",50,-7,7);
637 TH1F *hpz5 = new TH1F("hpz5","z pull - bin 5",50,-7,7);
638
639
640 TH1F *hmult = new TH1F("hmult","hmult",50,-0.5,49.5);
641
fee9ef0d 642 TFile *in = new TFile(inname.Data());
820252d7 643 TNtuple *nt = (TNtuple*)in->Get("ntVtxRes");
fee9ef0d 644 nt->SetBranchAddress("triggered",&triggered);
645 nt->SetBranchAddress("ntracks",&nESDTrks);
646 nt->SetBranchAddress("nitstracks",&nTrks);
820252d7 647 nt->SetBranchAddress("nitstracks5or6",&nTrks5or6);
820252d7 648 nt->SetBranchAddress("diffX",&diffX);
649 nt->SetBranchAddress("diffY",&diffY);
650 nt->SetBranchAddress("diffZ",&diffZ);
651 nt->SetBranchAddress("diffXerr",&diffXerr);
652 nt->SetBranchAddress("diffYerr",&diffYerr);
653 nt->SetBranchAddress("diffZerr",&diffZerr);
fee9ef0d 654 nt->SetBranchAddress("multiplicity",&ntracklets);
655 nt->SetBranchAddress("zMC",&zMC);
656 nt->SetBranchAddress("nc",&nc);
820252d7 657 Int_t entries = (Int_t)nt->GetEntries();
658 Int_t nbytes=0;
659
660 for(Int_t i=0; i<entries; i++) {
661 nbytes += nt->GetEvent(i);
fee9ef0d 662 // only triggered events
43c9dae1 663 // if(triggered<0.5) continue;
fee9ef0d 664 // consider only events with zMCmin<|zMC|<zMCmax
665 if(TMath::Abs(zMC)<zMCmin || TMath::Abs(zMC)>zMCmax) continue;
666 //
667 bin = GetBinTracklets(ntracklets);
668 zbin = GetBinZ(zMC);
669 mult[bin] += ntracklets;
670 hmult->Fill(ntracklets);
820252d7 671 totevts[bin]++;
fee9ef0d 672 hntrackletsVSzMC->Fill(zMC,ntracklets);
673 hnitstracksVSzMC->Fill(zMC,nTrks);
674 hnitstracks5or6VSzMC->Fill(zMC,nTrks5or6);
43c9dae1 675 if(TMath::Abs(diffX) < maxdiffx) { // vtx found - closer than maxdiffx micron
820252d7 676 nEvVtx++;
fee9ef0d 677 mult2[bin] += ntracklets;
820252d7 678 vtx[bin]++;
679 tottracks[bin] += nTrks;
680 evts[bin]++;
681 if(bin==0) hdx0->Fill(diffX);
682 if(bin==1) hdx1->Fill(diffX);
683 if(bin==2) hdx2->Fill(diffX);
684 if(bin==3) hdx3->Fill(diffX);
685 if(bin==4) hdx4->Fill(diffX);
686 if(bin==5) hdx5->Fill(diffX);
687 if(bin==0) hdy0->Fill(diffY);
688 if(bin==1) hdy1->Fill(diffY);
689 if(bin==2) hdy2->Fill(diffY);
690 if(bin==3) hdy3->Fill(diffY);
691 if(bin==4) hdy4->Fill(diffY);
692 if(bin==5) hdy5->Fill(diffY);
693 if(bin==0) hdz0->Fill(diffZ);
694 if(bin==1) hdz1->Fill(diffZ);
695 if(bin==2) hdz2->Fill(diffZ);
696 if(bin==3) hdz3->Fill(diffZ);
697 if(bin==4) hdz4->Fill(diffZ);
698 if(bin==5) hdz5->Fill(diffZ);
699 hdx->Fill(diffX);
700 hdy->Fill(diffY);
701 hdz->Fill(diffZ);
fee9ef0d 702 if(nc>=ncminforzshift) hdz_z[zbin].Fill(diffZ);
703 hzMCVtx->Fill(zMC);
704 hTrackletsVtx->Fill(ntracklets);
820252d7 705 hTracksVtx->Fill(nTrks);
706 hTracksVtx5or6->Fill(nTrks5or6);
820252d7 707 hPullsx->Fill(diffXerr);
708 hPullsy->Fill(diffYerr);
709 hPullsz->Fill(diffZerr);
710 if(bin==0) hpx0->Fill(diffXerr);
711 if(bin==1) hpx1->Fill(diffXerr);
712 if(bin==2) hpx2->Fill(diffXerr);
713 if(bin==3) hpx3->Fill(diffXerr);
714 if(bin==4) hpx4->Fill(diffXerr);
715 if(bin==5) hpx5->Fill(diffXerr);
716 if(bin==0) hpy0->Fill(diffYerr);
717 if(bin==1) hpy1->Fill(diffYerr);
718 if(bin==2) hpy2->Fill(diffYerr);
719 if(bin==3) hpy3->Fill(diffYerr);
720 if(bin==4) hpy4->Fill(diffYerr);
721 if(bin==5) hpy5->Fill(diffYerr);
722 if(bin==0) hpz0->Fill(diffZerr);
723 if(bin==1) hpz1->Fill(diffZerr);
724 if(bin==2) hpz2->Fill(diffZerr);
725 if(bin==3) hpz3->Fill(diffZerr);
726 if(bin==4) hpz4->Fill(diffZerr);
727 if(bin==5) hpz5->Fill(diffZerr);
728 } else {
729 nEvNoVtx++;
fee9ef0d 730 hzMCNoVtx->Fill(zMC);
731 hTrackletsNoVtx->Fill(ntracklets);
820252d7 732 hTracksNoVtx->Fill(nTrks);
733 hTracksNoVtx5or6->Fill(nTrks5or6);
820252d7 734 }
735 all[bin]++;
736 }
737
738
739 for(Int_t k=0;k<6;k++) {
740 mult2[k] /= evts[k];
741 mult[k] /= totevts[k];
742 tottracks[k] /= evts[k];
743 probVtx[k] = (Double_t)vtx[k]/all[k];
744 eprobVtx[k] = (Double_t)vtx[k]/all[k]/all[k]+(Double_t)vtx[k]*vtx[k]/all[k]/all[k]/all[k];
745 eprobVtx[k] = TMath::Sqrt(eprobVtx[k]);
746 }
747
748 // calculate spread in multiplicity
749 for(Int_t i=0; i<entries; i++) {
750 nbytes += nt->GetEvent(i);
fee9ef0d 751 bin = GetBinTracklets(ntracklets);
752 sigmamult[bin] += (ntracklets-mult[bin])*(ntracklets-mult[bin])/totevts[bin];
753 if(diffX < 90000.) sigmamult2[bin] += (ntracklets-mult2[bin])*(ntracklets-mult2[bin])/evts[bin];
820252d7 754 }
755
756 for(Int_t k=0;k<6;k++) {
757 sigmamult[k] = TMath::Sqrt(sigmamult[k]);
758 sigmamult2[k] = TMath::Sqrt(sigmamult2[k]);
759 }
760
761 // fraction of found vertices
762 Float_t frac = (Float_t)nEvVtx/(nEvVtx+nEvNoVtx);
763 printf(" Fraction of events with vertex: %f\n",frac);
764
765
766 // FIT FUNCTIONS
767 TF1 *g = new TF1("g","gaus");
768 g->SetLineWidth(3);
769 TF1 *line = new TF1("line","pol1");
770 line->SetLineWidth(3);
fee9ef0d 771 TF1 *func = new TF1("func",FitFunc,2,80,2);
820252d7 772 func->SetLineWidth(1);
773
774 Double_t x1,y1,sigmafit;
775 Char_t message[100];
776
777 //
778 // DRAW PLOTS
779 //
780 gStyle->SetOptFit(1111);
781
fee9ef0d 782 // tracks VS zMC
783 TCanvas *c0 = new TCanvas("c0","c0",0,0,1000,400);
784 c0->Divide(3,1);
785 c0->cd(1);
786 hntrackletsVSzMC->SetMinimum(0);
787 hntrackletsVSzMC->SetMaximum(14);
788 hntrackletsVSzMC->Draw("profs");
789 c0->cd(2);
790 hnitstracksVSzMC->SetMinimum(0);
791 hnitstracksVSzMC->SetMaximum(14);
792 hnitstracksVSzMC->Draw("profs");
793 c0->cd(3);
794 hnitstracks5or6VSzMC->SetMinimum(0);
795 hnitstracks5or6VSzMC->SetMaximum(14);
796 hnitstracks5or6VSzMC->Draw("profs");
797
798
820252d7 799 // Tracks in ITS for events w/ and w/o vertex
800 TCanvas *c1 = new TCanvas("c1","c1",0,0,1000,500);
801 c1->Divide(2,1);
802 c1->cd(1);
803 hTracksVtx->Draw();
804 c1->cd(2);
805 hTracksNoVtx->Draw();
806
fee9ef0d 807 // probability vs ntracklets
808 TCanvas *c1a = new TCanvas("c1a","c1a",0,0,500,500);
809 TH1F *hTotTracklets = (TH1F*)hTrackletsVtx->Clone("hTotTracklets");
810 hTotTracklets->Add(hTrackletsNoVtx);
811 TH1F *hProbTracklets = (TH1F*)hTrackletsVtx->Clone("hProbTracklets");
812 hProbTracklets->Divide(hTotTracklets);
813 hProbTracklets->SetYTitle("Probability");
814 hProbTracklets->SetTitle("Probability to find the vertex");
815 hProbTracklets->Draw();
816
820252d7 817 // probability vs ntracks
818 TCanvas *c1b = new TCanvas("c1b","c1b",0,0,500,500);
819 TH1F *hTotTracks = (TH1F*)hTracksVtx->Clone("hTotTracks");
820 hTotTracks->Add(hTracksNoVtx);
821 TH1F *hProbTracks = (TH1F*)hTracksVtx->Clone("hProbTracks");
822 hProbTracks->Divide(hTotTracks);
823 hProbTracks->SetYTitle("Probability");
824 hProbTracks->SetTitle("Probability to find the vertex");
825 hProbTracks->Draw();
826
827 TCanvas *c1c = new TCanvas("c1c","c1c",0,0,500,500);
828 TH1F *hTotTracks5or6 = (TH1F*)hTracksVtx5or6->Clone("hTotTracks5or6");
829 hTotTracks5or6->Add(hTracksNoVtx5or6);
830 TH1F *hProbTracks5or6 = (TH1F*)hTracksVtx5or6->Clone("hProbTracks5or6");
831 hProbTracks5or6->Divide(hTotTracks5or6);
832 hProbTracks5or6->SetYTitle("Probability");
833 hProbTracks5or6->SetTitle("Probability to find the vertex");
834 hProbTracks5or6->Draw();
835
fee9ef0d 836
837 // probability vs zMC
838 TCanvas *c1e = new TCanvas("c1e","c1e",0,0,500,500);
839 TH1F *hTotzMC = (TH1F*)hzMCVtx->Clone("hTotzMC");
840 hTotzMC->Add(hzMCNoVtx);
841 TH1F *hProbzMC = (TH1F*)hzMCVtx->Clone("hProbzMC");
842 hProbzMC->Divide(hTotzMC);
843 hProbzMC->SetYTitle("Probability");
844 hProbzMC->SetTitle("Probability to find the vertex");
845 hProbzMC->Draw();
820252d7 846
847 // Global resolutions
848 TCanvas *c2 = new TCanvas("c2","c2",0,0,1000,400);
849 c2->Divide(3,1);
850 c2->cd(1);
851 hdx->Draw();
852 g->SetRange(-1.*hdx->GetRMS(),+1.*hdx->GetRMS());
853 hdx->Fit("g","R");
854 c2->cd(2);
855 hdy->Draw();
856 g->SetRange(-1.*hdy->GetRMS(),+1.*hdy->GetRMS());
857 hdy->Fit("g","R");
858 c2->cd(3);
859 hdz->Draw();
860 g->SetRange(-1.*hdz->GetRMS(),+1.*hdz->GetRMS());
861 hdz->Fit("g","R");
862
863 // Pulls
864 TCanvas *c4 = new TCanvas("c4","c4",0,0,1000,400);
865 c4->Divide(3,1);
866 c4->cd(1);
867 hPullsx->Draw();
868 g->SetRange(-3.*hPullsx->GetRMS(),+3.*hPullsx->GetRMS());
869 hPullsx->Fit("g","R");
870 c4->cd(2);
871 hPullsy->Draw();
872 g->SetRange(-3.*hPullsy->GetRMS(),+3.*hPullsy->GetRMS());
873 hPullsy->Fit("g","R");
874 c4->cd(3);
875 hPullsz->Draw();
876 g->SetRange(-3.*hPullsz->GetRMS(),+3.*hPullsz->GetRMS());
877 hPullsz->Fit("g","R");
878
879 // probability VS multiplicity
880 TCanvas *c5 = new TCanvas("c5","c5",0,0,600,500);
fee9ef0d 881 TH2F *fr5 = new TH2F("fr5","",2,0,80,2,0,1.1);
882 fr5->SetXTitle("SPD tracklets");
820252d7 883 fr5->SetYTitle("Probability");
884 fr5->Draw();
885 TGraphErrors *gr5 = new TGraphErrors(6,mult,probVtx,sigmamult,eprobVtx);
886 gr5->Draw("p");
887 gr5->SetMarkerStyle(21);
888
889 //
890 // resolutions VS multiplicity
891 //
892 // X
893 TCanvas *ax = new TCanvas("ax","ax",0,0,900,700);
894 ax->Divide(3,2);
895 ax->cd(1);
896 hdx0->Draw();
897 g->SetRange(-3.*hdx0->GetRMS(),+3.*hdx0->GetRMS());
898 hdx0->Fit("g","Q");
899 xres[0]=g->GetParameter(2);
900 exres[0]=g->GetParError(2);
901 ax->cd(2);
902 hdx1->Draw();
903 g->SetRange(-3.*hdx1->GetRMS(),+3.*hdx1->GetRMS());
904 hdx1->Fit("g","Q");
905 xres[1]=g->GetParameter(2);
906 exres[1]=g->GetParError(2);
907 ax->cd(3);
908 hdx2->Draw();
909 g->SetRange(-3.*hdx2->GetRMS(),+3.*hdx2->GetRMS());
910 hdx2->Fit("g","Q");
911 xres[2]=g->GetParameter(2);
912 exres[2]=g->GetParError(2);
913 ax->cd(4);
914 hdx3->Draw();
915 g->SetRange(-3.*hdx3->GetRMS(),+3.*hdx3->GetRMS());
916 hdx3->Fit("g","Q");
917 xres[3]=g->GetParameter(2);
918 exres[3]=g->GetParError(2);
919 ax->cd(5);
920 hdx4->Draw();
921 g->SetRange(-3.*hdx4->GetRMS(),+3.*hdx4->GetRMS());
922 hdx4->Fit("g","Q");
923 xres[4]=g->GetParameter(2);
924 exres[4]=g->GetParError(2);
925 ax->cd(6);
926 hdx5->Draw();
927 g->SetRange(-3.*hdx5->GetRMS(),+3.*hdx5->GetRMS());
928 hdx5->Fit("g","Q");
929 xres[5]=g->GetParameter(2);
930 exres[5]=g->GetParError(2);
931 // Y
932 TCanvas *ay = new TCanvas("ay","ay",0,0,900,700);
933 ay->Divide(3,2);
934 ay->cd(1);
935 hdy0->Draw();
936 g->SetRange(-3.*hdy0->GetRMS(),+3.*hdy0->GetRMS());
937 hdy0->Fit("g","Q");
938 yres[0]=g->GetParameter(2);
939 eyres[0]=g->GetParError(2);
940 ay->cd(2);
941 hdy1->Draw();
942 g->SetRange(-3.*hdy1->GetRMS(),+3.*hdy1->GetRMS());
943 hdy1->Fit("g","Q");
944 yres[1]=g->GetParameter(2);
945 eyres[1]=g->GetParError(2);
946 ay->cd(3);
947 hdy2->Draw();
948 g->SetRange(-3.*hdy2->GetRMS(),+3.*hdy2->GetRMS());
949 hdy2->Fit("g","Q");
950 yres[2]=g->GetParameter(2);
951 eyres[2]=g->GetParError(2);
952 ay->cd(4);
953 hdy3->Draw();
954 g->SetRange(-3.*hdy3->GetRMS(),+3.*hdy3->GetRMS());
955 hdy3->Fit("g","Q");
956 yres[3]=g->GetParameter(2);
957 eyres[3]=g->GetParError(2);
958 ay->cd(5);
959 hdy4->Draw();
960 g->SetRange(-3.*hdy4->GetRMS(),+3.*hdy4->GetRMS());
961 hdy4->Fit("g","Q");
962 yres[4]=g->GetParameter(2);
963 eyres[4]=g->GetParError(2);
964 ay->cd(6);
965 hdy5->Draw();
966 g->SetRange(-3.*hdy5->GetRMS(),+3.*hdy5->GetRMS());
967 hdy5->Fit("g","Q");
968 yres[5]=g->GetParameter(2);
969 eyres[5]=g->GetParError(2);
970 // Z
971 TCanvas *az = new TCanvas("az","az",0,0,900,700);
972 az->Divide(3,2);
973 az->cd(1);
974 hdz0->Draw();
975 g->SetRange(-3.*hdz0->GetRMS(),+3.*hdz0->GetRMS());
976 hdz0->Fit("g","Q");
977 zres[0]=g->GetParameter(2);
978 ezres[0]=g->GetParError(2);
979 az->cd(2);
980 hdz1->Draw();
981 g->SetRange(-3.*hdz1->GetRMS(),+3.*hdz1->GetRMS());
982 hdz1->Fit("g","Q");
983 zres[1]=g->GetParameter(2);
984 ezres[1]=g->GetParError(2);
985 az->cd(3);
986 hdz2->Draw();
987 g->SetRange(-3.*hdz2->GetRMS(),+3.*hdz2->GetRMS());
988 hdz2->Fit("g","Q");
989 zres[2]=g->GetParameter(2);
990 ezres[2]=g->GetParError(2);
991 az->cd(4);
992 hdz3->Draw();
993 g->SetRange(-3.*hdz3->GetRMS(),+3.*hdz3->GetRMS());
994 hdz3->Fit("g","Q");
995 zres[3]=g->GetParameter(2);
996 ezres[3]=g->GetParError(2);
997 az->cd(5);
998 hdz4->Draw();
999 g->SetRange(-3.*hdz4->GetRMS(),+3.*hdz4->GetRMS());
1000 hdz4->Fit("g","Q");
1001 zres[4]=g->GetParameter(2);
1002 ezres[4]=g->GetParError(2);
1003 az->cd(6);
1004 hdz5->Draw();
1005 g->SetRange(-3.*hdz5->GetRMS(),+3.*hdz5->GetRMS());
1006 hdz5->Fit("g","Q");
1007 zres[5]=g->GetParameter(2);
1008 ezres[5]=g->GetParError(2);
1009
1010
1011 //
1012 // pulls VS multiplicity
1013 //
1014 // X
1015 TCanvas *bx = new TCanvas("bx","bx",0,0,900,700);
1016 bx->Divide(3,2);
1017 bx->cd(1);
1018 hpx0->Draw();
1019 g->SetRange(-3.*hpx0->GetRMS(),+3.*hpx0->GetRMS());
1020 hpx0->Fit("g","Q");
1021 xpull[0]=g->GetParameter(2);
1022 expull[0]=g->GetParError(2);
1023 bx->cd(2);
1024 hpx1->Draw();
1025 g->SetRange(-3.*hpx1->GetRMS(),+3.*hpx1->GetRMS());
1026 hpx1->Fit("g","Q");
1027 xpull[1]=g->GetParameter(2);
1028 expull[1]=g->GetParError(2);
1029 bx->cd(3);
1030 hpx2->Draw();
1031 g->SetRange(-3.*hpx2->GetRMS(),+3.*hpx2->GetRMS());
1032 hpx2->Fit("g","Q");
1033 xpull[2]=g->GetParameter(2);
1034 expull[2]=g->GetParError(2);
1035 bx->cd(4);
1036 hpx3->Draw();
1037 g->SetRange(-3.*hpx3->GetRMS(),+3.*hpx3->GetRMS());
1038 hpx3->Fit("g","Q");
1039 xpull[3]=g->GetParameter(2);
1040 expull[3]=g->GetParError(2);
1041 bx->cd(5);
1042 hpx4->Draw();
1043 g->SetRange(-3.*hpx4->GetRMS(),+3.*hpx4->GetRMS());
1044 hpx4->Fit("g","Q");
1045 xpull[4]=g->GetParameter(2);
1046 expull[4]=g->GetParError(2);
1047 bx->cd(6);
1048 hpx5->Draw();
1049 g->SetRange(-3.*hpx5->GetRMS(),+3.*hpx5->GetRMS());
1050 hpx5->Fit("g","Q");
1051 xpull[5]=g->GetParameter(2);
1052 expull[5]=g->GetParError(2);
1053 // Y
1054 TCanvas *by = new TCanvas("by","by",0,0,900,700);
1055 by->Divide(3,2);
1056 by->cd(1);
1057 hpy0->Draw();
1058 g->SetRange(-3.*hpy0->GetRMS(),+3.*hpy0->GetRMS());
1059 hpy0->Fit("g","Q");
1060 ypull[0]=g->GetParameter(2);
1061 eypull[0]=g->GetParError(2);
1062 by->cd(2);
1063 hpy1->Draw();
1064 g->SetRange(-3.*hpy1->GetRMS(),+3.*hpy1->GetRMS());
1065 hpy1->Fit("g","Q");
1066 ypull[1]=g->GetParameter(2);
1067 eypull[1]=g->GetParError(2);
1068 by->cd(3);
1069 hpy2->Draw();
1070 g->SetRange(-3.*hpy2->GetRMS(),+3.*hpy2->GetRMS());
1071 hpy2->Fit("g","Q");
1072 ypull[2]=g->GetParameter(2);
1073 eypull[2]=g->GetParError(2);
1074 by->cd(4);
1075 hpy3->Draw();
1076 g->SetRange(-3.*hpy3->GetRMS(),+3.*hpy3->GetRMS());
1077 hpy3->Fit("g","Q");
1078 ypull[3]=g->GetParameter(2);
1079 eypull[3]=g->GetParError(2);
1080 by->cd(5);
1081 hpy4->Draw();
1082 g->SetRange(-3.*hpy4->GetRMS(),+3.*hpy4->GetRMS());
1083 hpy4->Fit("g","Q");
1084 ypull[4]=g->GetParameter(2);
1085 eypull[4]=g->GetParError(2);
1086 by->cd(6);
1087 hpy5->Draw();
1088 g->SetRange(-3.*hpy5->GetRMS(),+3.*hpy5->GetRMS());
1089 hpy5->Fit("g","Q");
1090 ypull[5]=g->GetParameter(2);
1091 eypull[5]=g->GetParError(2);
1092 // Z
1093 TCanvas *bz = new TCanvas("bz","bz",0,0,900,700);
1094 bz->Divide(3,2);
1095 bz->cd(1);
1096 hpz0->Draw();
1097 g->SetRange(-3.*hpz0->GetRMS(),+3.*hpz0->GetRMS());
1098 hpz0->Fit("g","Q");
1099 zpull[0]=g->GetParameter(2);
1100 ezpull[0]=g->GetParError(2);
1101 bz->cd(2);
1102 hpz1->Draw();
1103 g->SetRange(-3.*hpz1->GetRMS(),+3.*hpz1->GetRMS());
1104 hpz1->Fit("g","Q");
1105 zpull[1]=g->GetParameter(2);
1106 ezpull[1]=g->GetParError(2);
1107 bz->cd(3);
1108 hpz2->Draw();
1109 g->SetRange(-3.*hpz2->GetRMS(),+3.*hpz2->GetRMS());
1110 hpz2->Fit("g","Q");
1111 zpull[2]=g->GetParameter(2);
1112 ezpull[2]=g->GetParError(2);
1113 bz->cd(4);
1114 hpz3->Draw();
1115 g->SetRange(-3.*hpz3->GetRMS(),+3.*hpz3->GetRMS());
1116 hpz3->Fit("g","Q");
1117 zpull[3]=g->GetParameter(2);
1118 ezpull[3]=g->GetParError(2);
1119 bz->cd(5);
1120 hpz4->Draw();
1121 g->SetRange(-3.*hpz4->GetRMS(),+3.*hpz4->GetRMS());
1122 hpz4->Fit("g","Q");
1123 zpull[4]=g->GetParameter(2);
1124 ezpull[4]=g->GetParError(2);
1125 bz->cd(6);
1126 hpz5->Draw();
1127 g->SetRange(-3.*hpz5->GetRMS(),+3.*hpz5->GetRMS());
1128 hpz5->Fit("g","Q");
1129 zpull[5]=g->GetParameter(2);
1130 ezpull[5]=g->GetParError(2);
1131
1132 gStyle->SetOptFit(0);
1133
fee9ef0d 1134 // res VS ntracklets
820252d7 1135 TCanvas *c6 = new TCanvas("c6","c6",0,0,600,500);
43c9dae1 1136 TH2F *fr6 = new TH2F("fr6","",2,0,80,2,0,240*rangefactor);
fee9ef0d 1137 fr6->SetXTitle("SPD tracklets");
820252d7 1138 fr6->SetYTitle("#sigma [#mu m]");
1139 fr6->Draw();
1140 sigmamult2[0]=sigmamult2[1]=sigmamult2[2]=sigmamult2[3]=sigmamult2[4]=sigmamult2[5]=0.;
1141 TGraphErrors *gr6x = new TGraphErrors(6,mult2,xres,sigmamult2,exres);
1142 gr6x->Draw("p");
1143 gr6x->SetMarkerStyle(22);
fee9ef0d 1144 //gr6x->Fit("func","E,R");
820252d7 1145 TGraphErrors *gr6z = new TGraphErrors(6,mult2,zres,sigmamult2,ezres);
1146 gr6z->Draw("p");
1147 gr6z->SetMarkerStyle(26);
fee9ef0d 1148 //gr6z->Fit("func","E,R");
820252d7 1149 TLegend *leg6 = new TLegend(0.6,0.75,0.9,0.9);
1150 leg6->AddEntry(gr6x,"x/y coordinate","p");
1151 leg6->AddEntry(gr6z,"z coordinate","p");
1152 leg6->SetFillStyle(0);
1153 leg6->SetBorderSize(0);
1154 leg6->Draw();
1155
fee9ef0d 1156 // pull VS ntracklets
820252d7 1157 TCanvas *c8 = new TCanvas("c8","c8",0,0,600,500);
fee9ef0d 1158 TH2F *fr8 = new TH2F("fr8","",2,0,80,2,0,2.);
1159 fr8->SetXTitle("SPD tracklets");
820252d7 1160 fr8->SetYTitle("pull");
1161 fr8->Draw();
1162 TGraphErrors *gr8x = new TGraphErrors(6,mult2,xpull,sigmamult2,expull);
1163 gr8x->Draw("p");
1164 gr8x->SetMarkerStyle(22);
1165 TGraphErrors *gr8z = new TGraphErrors(6,mult2,zpull,sigmamult2,ezpull);
1166 gr8z->Draw("p");
1167 gr8z->SetMarkerStyle(26);
1168 TLegend *leg8 = new TLegend(0.6,0.75,0.9,0.9);
1169 leg8->AddEntry(gr6x,"x/y coordinate","p");
1170 leg8->AddEntry(gr6z,"z coordinate","p");
1171 leg8->SetFillStyle(0);
1172 leg8->SetBorderSize(0);
1173 leg8->Draw();
fee9ef0d 1174 TLine *l8 = new TLine(0,1,80,1);
820252d7 1175 l8->SetLineStyle(2);
1176 l8->Draw();
1177
fee9ef0d 1178
1179 // mean z res VS zMC
1180 Float_t zmc[15]={-14,-12,-10,-8,-6,-4,-2,0,2,4,6,8,10,12,14};
1181 Float_t ezmc[15];
1182 Float_t dzmean[15],edzmean[15];
1183 Int_t dummy;
1184 TCanvas *zz = new TCanvas("zz","zz",0,0,1000,1000);
1185 zz->Divide(5,3);
1186 for(l=0;l<15;l++) {
1187 zz->cd(l+1);
1188 hdz_z[l].Draw();
1189 g->SetRange(-3.*hdz_z[l].GetRMS(),+3.*hdz_z[l].GetRMS());
1190 hdz_z[l].Fit("g","Q");
1191 dzmean[l]=g->GetParameter(1);
1192 edzmean[l]=g->GetParError(1);
1193 ezmc[l]=1.;
1194 }
1195 TCanvas *zzz = new TCanvas("zzz","zzz",0,0,600,500);
43c9dae1 1196 TH2F *fr9 = new TH2F("fr9","",2,-15,15,2,-50*rangefactor,50*rangefactor);
fee9ef0d 1197 fr9->SetXTitle("z_{true} [cm]");
1198 fr9->SetYTitle("<z_{rec} - z_{true}> [#mu m]");
1199 fr9->Draw();
1200 TGraphErrors *grz = new TGraphErrors(15,zmc,dzmean,ezmc,edzmean);
1201 grz->Draw("pz");
1202 grz->SetMarkerStyle(20);
1203
1204
820252d7 1205 return;
1206}
1207//----------------------------------------------------------------------------
fee9ef0d 1208void SaveFigures(TString name="dummy") {
1209
1210 TString namefull;
1211
1212 namefull = name.Data();
1213 namefull.Prepend("plots/residualsx_");
1214 namefull.Append(".gif");
1215 ax->Print(namefull.Data());
1216 namefull = name.Data();
1217 namefull.Prepend("plots/residualsx_");
1218 namefull.Append(".ps");
1219 ax->Print(namefull.Data());
1220 namefull = name.Data();
1221 namefull.Prepend("plots/residualsx_");
1222 namefull.Append(".eps");
1223 ax->Print(namefull.Data());
1224 namefull = name.Data();
1225 namefull.Prepend("plots/residualsx_");
1226 namefull.Append(".root");
1227 ax->Print(namefull.Data());
1228
1229 namefull = name.Data();
1230 namefull.Prepend("plots/pullsx_");
1231 namefull.Append(".gif");
1232 bx->Print(namefull.Data());
1233 namefull = name.Data();
1234 namefull.Prepend("plots/pullsx_");
1235 namefull.Append(".ps");
1236 bx->Print(namefull.Data());
1237 namefull = name.Data();
1238 namefull.Prepend("plots/pullsx_");
1239 namefull.Append(".eps");
1240 bx->Print(namefull.Data());
1241 namefull = name.Data();
1242 namefull.Prepend("plots/pullsx_");
1243 namefull.Append(".root");
1244 bx->Print(namefull.Data());
1245
1246 namefull = name.Data();
1247 namefull.Prepend("plots/residualsz_");
1248 namefull.Append(".gif");
1249 az->Print(namefull.Data());
1250 namefull = name.Data();
1251 namefull.Prepend("plots/residualsz_");
1252 namefull.Append(".ps");
1253 az->Print(namefull.Data());
1254 namefull = name.Data();
1255 namefull.Prepend("plots/residualsz_");
1256 namefull.Append(".eps");
1257 az->Print(namefull.Data());
1258 namefull = name.Data();
1259 namefull.Prepend("plots/residualsz_");
1260 namefull.Append(".root");
1261 az->Print(namefull.Data());
1262
1263 namefull = name.Data();
1264 namefull.Prepend("plots/pullsz_");
1265 namefull.Append(".gif");
1266 bz->Print(namefull.Data());
1267 namefull = name.Data();
1268 namefull.Prepend("plots/pullsz_");
1269 namefull.Append(".ps");
1270 bz->Print(namefull.Data());
1271 namefull = name.Data();
1272 namefull.Prepend("plots/pullsz_");
1273 namefull.Append(".eps");
1274 bz->Print(namefull.Data());
1275 namefull = name.Data();
1276 namefull.Prepend("plots/pullsz_");
1277 namefull.Append(".root");
1278 bz->Print(namefull.Data());
1279
1280 namefull = name.Data();
1281 namefull.Prepend("plots/effVSntracklets_");
1282 namefull.Append(".gif");
1283 c1a->Print(namefull.Data());
1284 namefull = name.Data();
1285 namefull.Prepend("plots/effVSntracklets_");
1286 namefull.Append(".ps");
1287 c1a->Print(namefull.Data());
1288 namefull = name.Data();
1289 namefull.Prepend("plots/effVSntracklets_");
1290 namefull.Append(".eps");
1291 c1a->Print(namefull.Data());
1292 namefull = name.Data();
1293 namefull.Prepend("plots/effVSntracklets_");
1294 namefull.Append(".root");
1295 c1a->Print(namefull.Data());
1296
1297
1298 namefull = name.Data();
1299 namefull.Prepend("plots/effVSnitstracks_");
1300 namefull.Append(".gif");
1301 c1b->Print(namefull.Data());
1302 namefull = name.Data();
1303 namefull.Prepend("plots/effVSnitstracks_");
1304 namefull.Append(".ps");
1305 c1b->Print(namefull.Data());
1306 namefull = name.Data();
1307 namefull.Prepend("plots/effVSnitstracks_");
1308 namefull.Append(".eps");
1309 c1b->Print(namefull.Data());
1310 namefull = name.Data();
1311 namefull.Prepend("plots/effVSnitstracks_");
1312 namefull.Append(".root");
1313 c1b->Print(namefull.Data());
1314
1315
1316 namefull = name.Data();
1317 namefull.Prepend("plots/effVSnitstracks5or6_");
1318 namefull.Append(".gif");
1319 c1c->Print(namefull.Data());
1320 namefull = name.Data();
1321 namefull.Prepend("plots/effVSnitstracks5or6_");
1322 namefull.Append(".ps");
1323 c1c->Print(namefull.Data());
1324 namefull = name.Data();
1325 namefull.Prepend("plots/effVSnitstracks5or6_");
1326 namefull.Append(".eps");
1327 c1c->Print(namefull.Data());
1328 namefull = name.Data();
1329 namefull.Prepend("plots/effVSnitstracks5or6_");
1330 namefull.Append(".root");
1331 c1c->Print(namefull.Data());
1332
1333
1334 namefull = name.Data();
1335 namefull.Prepend("plots/effVSzMC_");
1336 namefull.Append(".gif");
1337 c1e->Print(namefull.Data());
1338 namefull = name.Data();
1339 namefull.Prepend("plots/effVSzMC_");
1340 namefull.Append(".ps");
1341 c1e->Print(namefull.Data());
1342 namefull = name.Data();
1343 namefull.Prepend("plots/effVSzMC_");
1344 namefull.Append(".eps");
1345 c1e->Print(namefull.Data());
1346 namefull = name.Data();
1347 namefull.Prepend("plots/effVSzMC_");
1348 namefull.Append(".root");
1349 c1e->Print(namefull.Data());
1350
1351
1352
1353 namefull = name.Data();
1354 namefull.Prepend("plots/effVSntracklets2_");
1355 namefull.Append(".gif");
1356 c5->Print(namefull.Data());
1357 namefull = name.Data();
1358 namefull.Prepend("plots/effVSntracklets2_");
1359 namefull.Append(".ps");
1360 c5->Print(namefull.Data());
1361 namefull = name.Data();
1362 namefull.Prepend("plots/effVSntracklets2_");
1363 namefull.Append(".eps");
1364 c5->Print(namefull.Data());
1365 namefull = name.Data();
1366 namefull.Prepend("plots/effVSntracklets2_");
1367 namefull.Append(".root");
1368 c5->Print(namefull.Data());
1369
1370
1371 namefull = name.Data();
1372 namefull.Prepend("plots/sigmaVSntracklets_");
1373 namefull.Append(".gif");
1374 c6->Print(namefull.Data());
1375 namefull = name.Data();
1376 namefull.Prepend("plots/sigmaVSntracklets_");
1377 namefull.Append(".ps");
1378 c6->Print(namefull.Data());
1379 namefull = name.Data();
1380 namefull.Prepend("plots/sigmaVSntracklets_");
1381 namefull.Append(".eps");
1382 c6->Print(namefull.Data());
1383 namefull = name.Data();
1384 namefull.Prepend("plots/sigmaVSntracklets_");
1385 namefull.Append(".root");
1386 c6->Print(namefull.Data());
1387
1388
1389 namefull = name.Data();
1390 namefull.Prepend("plots/pullsVSntracklets_");
1391 namefull.Append(".gif");
1392 c8->Print(namefull.Data());
1393 namefull = name.Data();
1394 namefull.Prepend("plots/pullsVSntracklets_");
1395 namefull.Append(".ps");
1396 c8->Print(namefull.Data());
1397 namefull = name.Data();
1398 namefull.Prepend("plots/pullsVSntracklets_");
1399 namefull.Append(".eps");
1400 c8->Print(namefull.Data());
1401 namefull = name.Data();
1402 namefull.Prepend("plots/pullsVSntracklets_");
1403 namefull.Append(".root");
1404 c8->Print(namefull.Data());
1405
1406
1407 return;
1408}
1409//----------------------------------------------------------------------------
1410void TestRmTracks(Int_t iev=0) {
1411
1412 if (gAlice) {
33c3c91a 1413 delete AliRunLoader::Instance();
fee9ef0d 1414 delete gAlice;
1415 gAlice=0;
1416 }
1417
1418 AliRunLoader *rl = AliRunLoader::Open("galice.root");
1419 if (rl == 0x0) {
1420 cerr<<"Can not open session"<<endl;
1421 return;
1422 }
1423 Int_t retval = rl->LoadgAlice();
1424 if (retval) {
1425 cerr<<"LoadgAlice returned error"<<endl;
1426 delete rl;
1427 return;
1428 }
1429 retval = rl->LoadHeader();
1430 if (retval) {
1431 cerr<<"LoadHeader returned error"<<endl;
1432 delete rl;
1433 return;
1434 }
1435 gAlice=rl->GetAliRun();
1436 // Get field from galice.root
1437 AliMagF *fiel = (AliMagF*)gAlice->Field();
1438 // Set the conversion constant between curvature and Pt
1439 AliTracker::SetFieldMap(fiel,kTRUE);
1440
1441 Double_t truevtx[3];
1442 Double_t vtx[3];
1443 Double_t errvtx[3];
1444 Double_t diffX,diffY,diffZ;
1445 Double_t diffXerr,diffYerr,diffZerr;
1446 Double_t multiplicity;
1447 Int_t nc;
1448
1449 // ----------- Create vertexer --------------------------
1450 AliESDVertex *initVertex = 0;
1451 TFile *invtx = new TFile("AliESDVertexMean.root");
1452 initVertex = (AliESDVertex*)invtx->Get("vtxmean");
1453 invtx->Close();
1454 delete invtx;
1455
28e4388a 1456 AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
fee9ef0d 1457 vertexer->SetVtxStart(initVertex);
1458 vertexer->SetOnlyFitter();
1459 vertexer->SetDebug(0); // set to 1 to see what it does
1460
1461 Float_t diamondxy[2];
1462 diamondxy[0]=initVertex->GetXv();
1463 diamondxy[1]=initVertex->GetYv();
1464 //----------------------------------------------------------
1465
1466 if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
1467 printf("AliESDs.root not found!\n");
1468 return;
1469 }
1470 TFile *fin = TFile::Open("AliESDs.root");
64cf164e 1471 AliESDEvent *esd = new AliESDEvent();
fee9ef0d 1472 TTree *esdTree = (TTree*)fin->Get("esdTree");
1473 if(!esdTree) return;
64cf164e 1474 esd->ReadFromTree(esdTree);
fee9ef0d 1475
1476 TArrayF o;
1477
1478 nc=0;
1479
1480
1481 esdTree->GetEvent(iev);
1482 Int_t ntracks = esd->GetNumberOfTracks();
1483 printf("Number of tracks in ESD %d : %d\n",iev,ntracks);
1484
1485 // example: do vertex without tracks 1 and 2 (AliESDtrack::GetID())
1486 //
1487 // 1: tell vertexer to skip those two tracks
1488 //
1489 Int_t nskip=2;
1490 Int_t skipped[2];
1491 skipped[0]=1;
1492 skipped[1]=2;
1493 vertexer->SetSkipTracks(nskip,skipped);
1494 AliESDVertex *vertex1 = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
1495 vertex1->PrintStatus();
1496 vertex1->PrintIndices();
1497 delete vertex1;
1498 vertex1 = 0;
1499 //
1500 // 2: do vertex with all tracks
1501 //
1502 skipped[0]=-99;
1503 skipped[1]=-99;
1504 vertexer->SetSkipTracks(nskip,skipped);
1505 AliESDVertex *vertex2 = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
1506 vertex2->PrintStatus();
1507 vertex2->PrintIndices();
1508 //
1509 // 3: now remove those two tracks from vertex2
1510 //
1511 TTree *trksTree = new TTree("trksTree","tree with tracks to be removed");
1512 AliESDtrack *esdTrack = 0;
1513 trksTree->Branch("tracks","AliESDtrack",&esdTrack);
1514 AliESDtrack *et = esd->GetTrack(1);
1515 esdTrack = new AliESDtrack(*et);
1516 trksTree->Fill();
1517 delete esdTrack;
1518 AliESDtrack *et = esd->GetTrack(2);
1519 esdTrack = new AliESDtrack(*et);
1520 trksTree->Fill();
1521 delete esdTrack;
28e4388a 1522 AliVertexerTracks vt(AliTracker::GetBz());
fee9ef0d 1523 AliESDVertex *vertex3 = vt.RemoveTracksFromVertex(vertex2,trksTree,diamondxy);
1524 vertex3->PrintStatus();
1525 vertex3->PrintIndices();
1526 delete vertex2;
1527 vertex2 = 0;
1528 delete vertex3;
1529 vertex3 = 0;
1530
1531
1532 delete esdTree;
1533 delete fin;
1534 vertexer =0;
1535 delete vertexer;
1536 delete initVertex;
1537
1538 return;
1539}
1540//----------------------------------------------------------------------------
820252d7 1541void AliComputeVtxMeanFromESD(TString file="AliESDs.root",
1542 Int_t mincontr=20) {
1543 //-----------------------------------------------------------------------
1544 // Compute weighted mean and cov. matrix from a set of AliESD
1545 //-----------------------------------------------------------------------
1546
1547 gStyle->SetOptStat(0);
1548 //gStyle->SetOptFit(0);
1549
1550
1551 Double_t vtx[3],covvtx[6];
1552 TTree *esdTree = 0;
64cf164e 1553 AliESDEvent *esd = new AliESDEvent();
820252d7 1554 AliESDVertex *vertex = 0;
1555 TString vtitle;
1556 Int_t nc,events,total=0;
1557
1558 Double_t avx=0.;
1559 Double_t avy=0.;
1560 Double_t avz=0.;
1561 Double_t wgtavx=0.;
1562 Double_t wgtavy=0.;
1563 Double_t wgtavz=0.;
1564 Double_t sum=0.;
1565 Double_t sumwgtx=0.;
1566 Double_t sumwgty=0.;
1567 Double_t sumwgtz=0.;
1568 Double_t rmsx=0;
1569 Double_t rmsy=0;
1570 Double_t rmsz=0;
1571 Double_t varx=0;
1572 Double_t vary=0;
1573 Double_t varz=0;
1574 Double_t covxy=0;
1575 Double_t covxz=0;
1576 Double_t covyz=0;
1577 Double_t eavx,eavy,eavz,ewgtavx,ewgtavy,ewgtavz;
1578
1579 TH1F* hx = new TH1F("hx","",200,-1000,1000);
1580 hx->SetXTitle("vertex x [#mu m]");
1581 hx->SetYTitle("events");
1582 TH1F* hy = new TH1F("hy","",200,-1000,1000);
1583 hy->SetXTitle("vertex y [#mu m]");
1584 hy->SetYTitle("events");
1585 TH1F* hz = new TH1F("hz","",200,-10,10);
1586 hz->SetXTitle("vertex z [cm]");
1587 hz->SetYTitle("events");
1588
1589
1590 for(Int_t i=1; i<100; i++) {
1591 TString inname="event.";
1592 inname+=i;
1593 inname.Append("/");
1594 inname.Append(file.Data());
1595 if(gSystem->AccessPathName(inname.Data())) continue;
1596 TFile *fin = TFile::Open(inname.Data());
1597 esdTree = (TTree*)fin->Get("esdTree");
1598 if(!esdTree) continue;
64cf164e 1599 esd->ReadFromTree(esdTree);
820252d7 1600 events = esdTree->GetEntries();
1601 printf("Number of events in ESD tree: %d\n",events);
1602 total += events;
1603
1604 for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
1605 esdTree->GetEvent(iev);
1606 vertex = (AliESDVertex*)esd->GetPrimaryVertex();
1607 if(!vertex) continue;
1608 nc = vertex->GetNContributors();
1609 if(nc<mincontr) continue;
1610 vertex->GetXYZ(vtx);
1611 vertex->GetSigmaXYZ(errvtx);
1612 //printf("%f %f %f %f %f %f\n",vtx[0],vtx[1],vtx[2],errvtx[0],errvtx[1],errvtx[2]);
1613
1614 avx += vtx[0];
1615 avy += vtx[1];
1616 avz += vtx[2];
1617 sum += 1.;
1618 wgtavx += vtx[0]/errvtx[0]/errvtx[0];
1619 wgtavy += vtx[1]/errvtx[1]/errvtx[1];
1620 wgtavz += vtx[2]/errvtx[2]/errvtx[2];
1621 sumwgtx += 1./errvtx[0]/errvtx[0];
1622 sumwgty += 1./errvtx[1]/errvtx[1];
1623 sumwgtz += 1./errvtx[2]/errvtx[2];
1624
1625 hx->Fill(1e4*vtx[0]);
1626 hy->Fill(1e4*vtx[1]);
1627 hz->Fill(vtx[2]);
1628
1629 vertex = 0;
1630 }
1631 esdTree=0;
1632 fin->Close();
1633 delete fin;
1634 }
1635
1636 avx /= sum;
1637 avy /= sum;
1638 avz /= sum;
1639 wgtavx /= sumwgtx;
1640 wgtavy /= sumwgty;
1641 wgtavz /= sumwgtz;
1642 ewgtavx = 1./TMath::Sqrt(sumwgtx);
1643 ewgtavy = 1./TMath::Sqrt(sumwgty);
1644 ewgtavz = 1./TMath::Sqrt(sumwgtz);
1645
1646
1647 for(Int_t i=1; i<100; i++) {
1648 TString inname="event.";
1649 inname+=i;
1650 inname.Append("/");
1651 inname.Append(file.Data());
1652 if(gSystem->AccessPathName(inname.Data())) continue;
1653 TFile *fin = TFile::Open(inname.Data());
1654 esdTree = (TTree*)fin->Get("esdTree");
1655 if(!esdTree) continue;
64cf164e 1656 esd->ReadFromTree(esdTree);
820252d7 1657 events = esdTree->GetEntries();
1658 for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
1659 esdTree->GetEvent(iev);
1660 vertex = (AliESDVertex*)esd->GetPrimaryVertex();
1661 if(!vertex) continue;
1662 nc = vertex->GetNContributors();
1663 if(nc<mincontr) continue;
1664 vertex->GetXYZ(vtx);
1665 vertex->GetSigmaXYZ(errvtx);
1666
1667 varx += (vtx[0]-avx)*(vtx[0]-avx);
1668 vary += (vtx[1]-avy)*(vtx[1]-avy);
1669 varz += (vtx[2]-avz)*(vtx[2]-avz);
1670 covxy += (vtx[0]-avx)*(vtx[1]-avy);
1671 covxz += (vtx[0]-avx)*(vtx[2]-avz);
1672 covyz += (vtx[1]-avy)*(vtx[2]-avz);
1673
1674 vertex = 0;
1675 }
1676 }
1677
1678 varx /= sum;
1679 vary /= sum;
1680 varz /= sum;
1681 covxy /= sum;
1682 covxz /= sum;
1683 covyz /= sum;
1684 rmsx = TMath::Sqrt(varx);
1685 rmsy = TMath::Sqrt(vary);
1686 rmsz = TMath::Sqrt(varz);
1687 eavx = rmsx/TMath::Sqrt(sum);
1688 eavy = rmsy/TMath::Sqrt(sum);
1689 eavz = rmsz/TMath::Sqrt(sum);
1690
1691
1692 printf("\n\nNumber of events: Total %d, Used %d\n",total,sum);
1693 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);
1694 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);
1695 printf("RMS:\n x = %f mum\n y = %f mum\n z = %f mum\n",1e4*rmsx,1e4*rmsy,1e4*rmsz);
1696
1697 TCanvas *c = new TCanvas("c","c",0,0,1000,500);
1698 c->Divide(3,1);
1699 c->cd(1);
1700 hx->Draw();
1701 TF1 *gx = new TF1("gx","gaus",-1000,1000);
1702 gx->SetLineColor(2);
1703 hx->Fit(gx);
1704 TF1 *gxx = (TF1*)gx->Clone("gxx");
1705 gxx->FixParameter(2,50.);
1706 gxx->SetLineStyle(2);
1707 gxx->Draw("same");
1708 c->cd(2);
1709 hy->Draw();
1710 TF1 *gy = new TF1("gy","gaus",-1000,1000);
1711 gy->SetLineColor(2);
1712 hy->Fit(gy);
1713 TF1 *gyy = (TF1*)gy->Clone("gyy");
1714 gyy->FixParameter(2,50.);
1715 gyy->SetLineStyle(2);
1716 gyy->Draw("same");
1717 c->cd(3);
1718 hz->Draw();
1719 TF1 *gz = new TF1("gz","gaus",-10,10);
1720 gz->SetLineColor(2);
1721 hz->Fit(gz);
1722 TF1 *gzz = (TF1*)gz->Clone("gzz");
1723 gzz->FixParameter(2,5.3);
1724 gzz->SetLineStyle(2);
1725 gzz->Draw("same");
1726
1727 vtx[0] = wgtavx;
1728 vtx[1] = wgtavy;
1729 vtx[2] = wgtavz;
1730 covvtx[0] = varx;
1731 covvtx[1] = covxy;
1732 covvtx[2] = vary;
1733 covvtx[3] = covxz;
1734 covvtx[4] = covyz;
1735 covvtx[5] = 5.3;
1736 AliESDVertex *vtxmean = new AliESDVertex(vtx,covvtx,0.,(Int_t)sum);
1737 vtxmean->SetTitle("Vertex from weighted average");
1738 TFile *out = new TFile("AliESDVertexMean.root","recreate");
1739 vtxmean->Write("vtxmean");
1740 out->Close();
1741
1742 return;
1743}
1744//----------------------------------------------------------------------------
1745void AliComputeVtxMeanFromTree(TString file="AliESDVertices.root",
1746 Int_t mincontr=20) {
1747 //-----------------------------------------------------------------------
1748 // Compute weighted mean and cov. matrix from a set of AliESDVertices
1749 //-----------------------------------------------------------------------
1750 gStyle->SetOptStat(0);
1751 //gStyle->SetOptFit(0);
1752
1753
1754 Double_t vtx[3],covvtx[6];
1755 TString vtitle;
1756 Int_t nc,total=0;
1757
1758 Double_t avx=0.;
1759 Double_t avy=0.;
1760 Double_t avz=0.;
1761 Double_t wgtavx=0.;
1762 Double_t wgtavy=0.;
1763 Double_t wgtavz=0.;
1764 Double_t sum=0.;
1765 Double_t sumwgtx=0.;
1766 Double_t sumwgty=0.;
1767 Double_t sumwgtz=0.;
1768 Double_t rmsx=0;
1769 Double_t rmsy=0;
1770 Double_t rmsz=0;
1771 Double_t eavx,eavy,eavz,ewgtavx,ewgtavy,ewgtavz;
1772
1773 TH1F* hx = new TH1F("hx","",200,-1000,1000);
1774 hx->SetXTitle("vertex x [#mu m]");
1775 hx->SetYTitle("events");
1776 TH1F* hy = new TH1F("hy","",200,-1000,1000);
1777 hy->SetXTitle("vertex y [#mu m]");
1778 hy->SetYTitle("events");
1779 TH1F* hz = new TH1F("hz","",200,-10,10);
1780 hz->SetXTitle("vertex z [cm]");
1781 hz->SetYTitle("events");
1782
1783
1784 if(gSystem->AccessPathName(file.Data())) continue;
1785 TFile *fin = TFile::Open(file.Data());
1786 TTree *tree = (TTree*)fin->Get("TreeV");
1787 AliESDVertex *vertex = 0;
1788 tree->SetBranchAddress("AliESDVertex",&vertex);
1789 Int_t entries = (Int_t)tree->GetEntries();
1790
1791 for(Int_t i=0; i<entries; i++) {
1792 total += 1;
1793 tree->GetEvent(i);
1794 nc = vertex->GetNContributors();
1795 if(nc<mincontr) continue;
1796 vertex->GetXYZ(vtx);
1797 vertex->GetSigmaXYZ(errvtx);
1798 //printf("%f %f %f %f %f %f\n",vtx[0],vtx[1],vtx[2],errvtx[0],errvtx[1],errvtx[2]);
1799
1800 avx += vtx[0];
1801 avy += vtx[1];
1802 avz += vtx[2];
1803 sum += 1.;
1804 wgtavx += vtx[0]/errvtx[0]/errvtx[0];
1805 wgtavy += vtx[1]/errvtx[1]/errvtx[1];
1806 wgtavz += vtx[2]/errvtx[2]/errvtx[2];
1807 sumwgtx += 1./errvtx[0]/errvtx[0];
1808 sumwgty += 1./errvtx[1]/errvtx[1];
1809 sumwgtz += 1./errvtx[2]/errvtx[2];
1810
1811 hx->Fill(1e4*vtx[0]);
1812 hy->Fill(1e4*vtx[1]);
1813 hz->Fill(vtx[2]);
1814
1815 vertex = 0;
1816 }
1817
1818 avx /= sum;
1819 avy /= sum;
1820 avz /= sum;
1821 wgtavx /= sumwgtx;
1822 wgtavy /= sumwgty;
1823 wgtavz /= sumwgtz;
1824 ewgtavx = 1./TMath::Sqrt(sumwgtx);
1825 ewgtavy = 1./TMath::Sqrt(sumwgty);
1826 ewgtavz = 1./TMath::Sqrt(sumwgtz);
1827
1828
1829 for(Int_t i=0; i<entries; i++) {
1830 tree->GetEvent(i);
1831 nc = vertex->GetNContributors();
1832 if(nc<mincontr) continue;
1833 vertex->GetXYZ(vtx);
1834 vertex->GetSigmaXYZ(errvtx);
1835
1836 varx += (vtx[0]-avx)*(vtx[0]-avx);
1837 vary += (vtx[1]-avy)*(vtx[1]-avy);
1838 varz += (vtx[2]-avz)*(vtx[2]-avz);
1839 covxy += (vtx[0]-avx)*(vtx[1]-avy);
1840 covxz += (vtx[0]-avx)*(vtx[2]-avz);
1841 covyz += (vtx[1]-avy)*(vtx[2]-avz);
1842
1843 vertex = 0;
1844 }
1845
1846 varx /= sum;
1847 vary /= sum;
1848 varz /= sum;
1849 covxy /= sum;
1850 covxz /= sum;
1851 covyz /= sum;
1852 rmsx = TMath::Sqrt(varx);
1853 rmsy = TMath::Sqrt(vary);
1854 rmsz = TMath::Sqrt(varz);
1855 eavx = rmsx/TMath::Sqrt(sum);
1856 eavy = rmsy/TMath::Sqrt(sum);
1857 eavz = rmsz/TMath::Sqrt(sum);
1858
1859
1860 printf("\n\nNumber of events: Total %d, Used %d\n",total,sum);
1861 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);
1862 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);
1863 printf("RMS:\n x = %f mum\n y = %f mum\n z = %f mum\n",1e4*rmsx,1e4*rmsy,1e4*rmsz);
1864
1865 TCanvas *c = new TCanvas("c","c",0,0,1000,500);
1866 c->Divide(3,1);
1867 c->cd(1);
1868 hx->Draw();
1869 TF1 *gx = new TF1("gx","gaus",-1000,1000);
1870 gx->SetLineColor(2);
1871 hx->Fit(gx);
1872 TF1 *gxx = (TF1*)gx->Clone("gxx");
1873 gxx->FixParameter(2,50.);
1874 gxx->SetLineStyle(2);
1875 gxx->Draw("same");
1876 c->cd(2);
1877 hy->Draw();
1878 TF1 *gy = new TF1("gy","gaus",-1000,1000);
1879 gy->SetLineColor(2);
1880 hy->Fit(gy);
1881 TF1 *gyy = (TF1*)gy->Clone("gyy");
1882 gyy->FixParameter(2,50.);
1883 gyy->SetLineStyle(2);
1884 gyy->Draw("same");
1885 c->cd(3);
1886 hz->Draw();
1887 TF1 *gz = new TF1("gz","gaus",-10,10);
1888 gz->SetLineColor(2);
1889 hz->Fit(gz);
1890 TF1 *gzz = (TF1*)gz->Clone("gzz");
1891 gzz->FixParameter(2,5.3);
1892 gzz->SetLineStyle(2);
1893 gzz->Draw("same");
1894
1895 vtx[0] = wgtavx;
1896 vtx[1] = wgtavy;
1897 vtx[2] = wgtavz;
1898 covvtx[0] = varx;
1899 covvtx[1] = covxy;
1900 covvtx[2] = vary;
1901 covvtx[3] = covxz;
1902 covvtx[4] = covyz;
1903 covvtx[5] = 5.3;
1904 AliESDVertex *vtxmean = new AliESDVertex(vtx,covvtx,0.,(Int_t)sum);
1905 vtxmean->SetTitle("Vertex from weighted average");
1906 TFile *out = new TFile("AliESDVertexMean.root","recreate");
1907 vtxmean->Write("vtxmean");
1908 out->Close();
1909
1910 return;
1911}
1912//---------------------------------------------------------------------------
43c9dae1 1913