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