STEER/AliVertexerTracksTest.C

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