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->GetXv();
1463   diamondxy[1]=initVertex->GetYv();
1464   //----------------------------------------------------------
1465
1466   if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
1467     printf("AliESDs.root not found!\n");
1468     return;
1469   }
1470   TFile *fin = TFile::Open("AliESDs.root");
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