Replaced Error(...); exit(1) by AliFatal in the AliCheb... classes
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracksTest.C
1 void AliVertexerTracksTest(TString outname="AliVertexerTracksTest.root",
2                            Bool_t tpconly=kTRUE,
3                            Bool_t useMeanVtx=kFALSE,
4                            Double_t nSigma=3.,
5                            Int_t itsrefit=1,
6                            Double_t maxd0z0=0.5,
7                            Int_t minitscls=5,
8                            Int_t mintracks=1,
9                            Bool_t onlyfit=kFALSE) {
10 //------------------------------------------------------------------------
11 // Run vertex reconstruction and store results in histos and ntuple
12 //------------------------------------------------------------------------
13 // WITHOUT Kinematics.root
14 //
15   AliGeomManager::LoadGeometry("geometry.root");
16
17   if (gAlice) {
18     delete AliRunLoader::Instance();
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();
41
42   //rl->LoadKinematics();
43
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;
59   Float_t multiplicity;
60   Float_t chi2red;
61   Int_t ntracklets;
62   Int_t nc;
63   Int_t triggered;
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
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");
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   }
92   AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
93   vertexer->SetDebug(1); // set to 1 to see what it does
94   vertexer->SetVtxStart(initVertex);
95   if(!useMeanVtx) vertexer->SetConstraintOff();
96   if(onlyfit) vertexer->SetOnlyFitter();
97   vertexer->SetNSigmad0(nSigma);
98   if(!itsrefit || tpconly) vertexer->SetITSrefitNotRequired();
99   vertexer->SetMinTracks(mintracks);
100   vertexer->SetMinITSClusters(minitscls);
101   vertexer->SetMaxd0z0(maxd0z0);
102   //cout<<" Nsigma:  "<<vertexer->GetNSigmad0()<<endl;
103   //----------------------------------------------------------
104  
105
106   if(gSystem->AccessPathName("AliESDs_itstpc.root",kFileExists)) {
107     printf("AliESDs_itstpc.root not found!\n"); 
108     return;
109   }
110   TFile *fin = TFile::Open("AliESDs_itstpc.root");
111   AliESDEvent *esd = new AliESDEvent();
112   TTree *esdTree = (TTree*)fin->Get("esdTree");
113   if(!esdTree) return;
114   esd->ReadFromTree(esdTree);
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
123     diffX=999999;
124     diffY=999999;
125     diffZ=999999;
126     diffXerr=999999;
127     diffYerr=999999;
128     diffZerr=999999;
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;
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;
154
155
156     esdTree->GetEvent(iev);
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
174     Int_t ntracks = esd->GetNumberOfTracks();
175     printf("Number of tracks in ESD %d   :   %d\n",iev,ntracks);
176     if(ntracks==0) { 
177       ntVtxRes->Fill(triggered,ntracks,nitstracks,nitstracks5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,nc,truevtx[2],chi2red); 
178       continue; 
179     }
180     
181     // VERTEX    
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
211
212     nc = (Int_t)vertex->GetNContributors();
213     if(nc>=2) chi2red = vertex->GetChi2toNDF(); 
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++;
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);
239
240
241     if(nc>=1) {
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
258     ntVtxRes->Fill(triggered,ntracks,nitstracks,nitstracks5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,(Float_t)nc,truevtx[2],chi2red);   
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
281 //------------------------------------------------------------------------
282 void 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) {
290     delete AliRunLoader::Instance();
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   }
340   AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
341   //vertexer->SetITSrefitNotRequired();
342   vertexer->SetVtxStart(initVertex);
343   if(!useMeanVtx) vertexer->SetConstraintOff();
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");
357   AliESDEvent *esd = new AliESDEvent();
358   TTree *esdTree = (TTree*)fin->Get("esdTree");
359   if(!esdTree) return;
360   esd->ReadFromTree(esdTree);
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 //----------------------------------------------------------------------------
429 Double_t FitFunc(Double_t *x,Double_t *par) {
430   return par[0]+par[1]/TMath::Sqrt(x[0]);
431 }
432 Int_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 }
440 Int_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 }
448 Int_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 //----------------------------------------------------------------------------
466 void PlotVtxRes(TString inname="AliVertexerTracksTest.root",
467                 Bool_t tpconly=kFALSE) {
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
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
486   Int_t   nEvVtx=0;
487   Int_t   nEvNoVtx=0;
488   Int_t   ev,nUsedTrks;
489   Float_t nESDTrks,nTrks,nTrks5or6,ntracklets;
490   Float_t diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,zMC,nc;
491   Float_t triggered;
492
493
494   //
495   // HISTOGRAMS
496   //
497   TH1F *hdx = new TH1F("hdx","",50,-600*rangefactor,600*rangefactor);
498   hdx->SetXTitle("#Delta x [#mu m]");
499   hdx->SetYTitle("events");
500   //
501   TH1F *hdy = new TH1F("hdy","",50,-600*rangefactor,600*rangefactor);
502   hdy->SetXTitle("#Delta y [#mu m]");
503   hdy->SetYTitle("events");
504   //
505   TH1F *hdz = new TH1F("hdz","",200,-600*rangefactor,600*rangefactor);
506   hdz->SetXTitle("#Delta z [#mu m]");
507   hdz->SetYTitle("events");
508   //
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   //
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);
536   hTracksVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(cls>=5)");
537   hTracksVtx5or6->SetYTitle("events");
538   //
539   TH1F *hTracksNoVtx5or6 = new TH1F("hTracksNoVtx5or6","events w/o vertex",51,-0.5,50.5);
540   hTracksNoVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(cls>=5)");
541   hTracksNoVtx5or6->SetYTitle("events");
542   //
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   //
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)>");
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};
593   Int_t    bin,zbin;
594   //
595   // OTHER HISTOGRAMS
596   //
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);
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
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
642   TFile *in = new TFile(inname.Data());
643   TNtuple *nt = (TNtuple*)in->Get("ntVtxRes");
644   nt->SetBranchAddress("triggered",&triggered);
645   nt->SetBranchAddress("ntracks",&nESDTrks);
646   nt->SetBranchAddress("nitstracks",&nTrks);
647   nt->SetBranchAddress("nitstracks5or6",&nTrks5or6);
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);
654   nt->SetBranchAddress("multiplicity",&ntracklets);
655   nt->SetBranchAddress("zMC",&zMC);
656   nt->SetBranchAddress("nc",&nc);
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);
662     // only triggered events
663     //    if(triggered<0.5) continue;
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);
671     totevts[bin]++;
672     hntrackletsVSzMC->Fill(zMC,ntracklets);
673     hnitstracksVSzMC->Fill(zMC,nTrks);
674     hnitstracks5or6VSzMC->Fill(zMC,nTrks5or6);
675     if(TMath::Abs(diffX) < maxdiffx) { // vtx found - closer than maxdiffx micron
676       nEvVtx++;
677       mult2[bin] += ntracklets;
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);
702       if(nc>=ncminforzshift) hdz_z[zbin].Fill(diffZ);
703       hzMCVtx->Fill(zMC);
704       hTrackletsVtx->Fill(ntracklets);
705       hTracksVtx->Fill(nTrks);
706       hTracksVtx5or6->Fill(nTrks5or6);
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++;
730       hzMCNoVtx->Fill(zMC);
731       hTrackletsNoVtx->Fill(ntracklets);
732       hTracksNoVtx->Fill(nTrks);
733       hTracksNoVtx5or6->Fill(nTrks5or6);
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);
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];
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);
771   TF1 *func = new TF1("func",FitFunc,2,80,2);
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   
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
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
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
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
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();
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);
881   TH2F *fr5 = new TH2F("fr5","",2,0,80,2,0,1.1); 
882   fr5->SetXTitle("SPD tracklets");
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
1134   // res VS ntracklets
1135   TCanvas *c6 = new TCanvas("c6","c6",0,0,600,500);
1136   TH2F *fr6 = new TH2F("fr6","",2,0,80,2,0,240*rangefactor); 
1137   fr6->SetXTitle("SPD tracklets");
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);
1144   //gr6x->Fit("func","E,R");
1145   TGraphErrors *gr6z = new TGraphErrors(6,mult2,zres,sigmamult2,ezres);
1146   gr6z->Draw("p");
1147   gr6z->SetMarkerStyle(26);
1148   //gr6z->Fit("func","E,R");
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
1156   // pull VS ntracklets
1157   TCanvas *c8 = new TCanvas("c8","c8",0,0,600,500);
1158   TH2F *fr8 = new TH2F("fr8","",2,0,80,2,0,2.); 
1159   fr8->SetXTitle("SPD tracklets");
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();
1174   TLine *l8 = new TLine(0,1,80,1);
1175   l8->SetLineStyle(2);
1176   l8->Draw();
1177
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);
1196   TH2F *fr9 = new TH2F("fr9","",2,-15,15,2,-50*rangefactor,50*rangefactor); 
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
1205   return;
1206 }
1207 //----------------------------------------------------------------------------
1208 void 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 //----------------------------------------------------------------------------
1410 void TestRmTracks(Int_t iev=0) {
1411
1412   if (gAlice) {
1413     delete AliRunLoader::Instance();
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
1456   AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
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->GetX();
1463   diamondxy[1]=initVertex->GetY();
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");
1471   AliESDEvent *esd = new AliESDEvent();
1472   TTree *esdTree = (TTree*)fin->Get("esdTree");
1473   if(!esdTree) return;
1474   esd->ReadFromTree(esdTree);
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;
1522   AliVertexerTracks vt(AliTracker::GetBz());
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 //----------------------------------------------------------------------------
1541 void 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;
1553   AliESDEvent *esd = new AliESDEvent();
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;
1599     esd->ReadFromTree(esdTree);
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;
1656     esd->ReadFromTree(esdTree);
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 //----------------------------------------------------------------------------
1745 void 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 //---------------------------------------------------------------------------
1913