STEER/AliVertexerTracksTest.C

   1 //------------------------------------------------------------------------
   2 // Test macro for AliVertexerTracks.
   3 // Contains several functions.
   4 //
   5 // A. Dainese - INFN Legnaro
   6 //------------------------------------------------------------------------
   7 void AliVertexerTracksTest(Double_t nSigma=3.,
   8                            Bool_t useMeanVtx=kTRUE,
   9                            TString outname="AliVertexerTracksTest.root") {
  10 //------------------------------------------------------------------------
  11 // Run vertex reconstruction and store results in histos and ntuple
  12 //------------------------------------------------------------------------
  13
  14   if (gAlice) {
  15     delete gAlice->GetRunLoader();
  16     delete gAlice;
  17     gAlice=0;
  18   }
  19
  20   AliRunLoader *rl = AliRunLoader::Open("galice.root");
  21   if (rl == 0x0) {
  22     cerr<<"Can not open session"<<endl;
  23     return;
  24   }
  25   Int_t retval = rl->LoadgAlice();
  26   if (retval) {
  27     cerr<<"LoadgAlice returned error"<<endl;
  28     delete rl;
  29     return;
  30   }
  31   retval = rl->LoadHeader();
  32   if (retval) {
  33     cerr<<"LoadHeader returned error"<<endl;
  34     delete rl;
  35     return;
  36   }
  37   gAlice=rl->GetAliRun();
  38   rl->LoadKinematics();
  39   // Get field from galice.root
  40   AliMagF *fiel = (AliMagF*)gAlice->Field();
  41   // Set the conversion constant between curvature and Pt
  42   AliTracker::SetFieldMap(fiel,kTRUE);
  43
  44   //**** Switch on the PID class
  45   AliPID pid;
  46
  47   Int_t nEvents = (Int_t)rl->GetNumberOfEvents();
  48
  49   Double_t truevtx[3];
  50   Double_t vtx[3];
  51   Double_t errvtx[3];
  52   Double_t diffX,diffY,diffZ;
  53   Double_t diffXerr,diffYerr,diffZerr;
  54   Double_t multiplicity;
  55   Int_t nc;
  56
  57   // Histograms
  58   TH1F *hdiffX = new TH1F("hdiffX","x_{found} - x_{true} [#mum]",1000,-5000,5000);
  59   TH1F *hdiffY = new TH1F("hdiffY","y_{found} - y_{true} [#mum]",1000,-5000,5000);
  60   TH1F *hdiffZ = new TH1F("hdiffZ","z_{found} - z_{true} [#mum]",1000,-5000,5000);
  61
  62   TH1F *hdiffXerr = new TH1F("hdiffXerr","#Delta x/#sigma_{x}",100,-5,5);
  63   TH1F *hdiffYerr = new TH1F("hdiffYerr","#Delta y/#sigma_{y}",100,-5,5);
  64   TH1F *hdiffZerr = new TH1F("hdiffZerr","#Delta z/#sigma_{z}",100,-5,5);
  65
  66   //TNtple
  67   TNtuple *ntVtxRes = new TNtuple("ntVtxRes","Vertex Resolution's Ntupla with multiplicity","ntracks:nitstracks5or6:nitstracksFromStrange:nitstracksFromStrange5or6:diffX:diffY:diffZ:diffXerr:diffYerr:diffZerr:multiplicity:nc");
  68
  69
  70
  71   // -----------   Create vertexer --------------------------
  72   AliESDVertex *initVertex = 0;
  73   if(useMeanVtx) {
  74     // open the mean vertex
  75     TFile *invtx = new TFile("AliESDVertexMean.root");
  76     initVertex = (AliESDVertex*)invtx->Get("vtxmean");
  77     invtx->Close();
  78     delete invtx;
  79   } else {
  80     Double_t pos[3]={0.5,0.5,0.};
  81     Double_t err[3]={3.,3.,5.};
  82     initVertex = new AliESDVertex(pos,err);
  83   }
  84   AliVertexerTracks *vertexer = new AliVertexerTracks();
  85   vertexer->SetITSrefitNotRequired();
  86   vertexer->SetVtxStart(initVertex);
  87   vertexer->SetMinTracks(2);
  88   vertexer->SetNSigmad0(nSigma);
  89   //cout<<" Nsigma:  "<<vertexer->GetNSigmad0()<<endl;
  90   vertexer->SetDebug(1); // set to 1 to see what it does
  91   //----------------------------------------------------------
  92
  93   if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
  94     printf("AliESDs.root not found!\n");
  95     return;
  96   }
  97   TFile *fin = TFile::Open("AliESDs.root");
  98   TTree *esdTree = (TTree*)fin->Get("esdTree");
  99   if(!esdTree) return;
 100   AliESD *esd = 0;
 101   esdTree->SetBranchAddress("ESD",&esd);
 102   Int_t events = esdTree->GetEntries();
 103   printf("Number of events in ESD tree: %d\n",events);
 104
 105   TArrayF o;
 106
 107   for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
 108     printf("-------------- EVENT   %d --------------------\n",iev);
 109
 110     diffX=99999;
 111     diffY=99999;
 112     diffZ=99999;
 113     diffXerr=99999;
 114     diffYerr=99999;
 115     diffZerr=99999;
 116     nc=0;
 117
 118     //=================================================
 119     //         LOOK IN THE SIMULATION
 120     // true vertex position
 121     Int_t npart = (Int_t)gAlice->GetEvent(iev);
 122     printf("particles  %d\n",npart);
 123     AliHeader *header = (AliHeader*)rl->GetHeader();
 124     AliGenPythiaEventHeader *pyheader = (AliGenPythiaEventHeader*)header->GenEventHeader();
 125     pyheader->PrimaryVertex(o);
 126     truevtx[0] = o[0];
 127     truevtx[1] = o[1];
 128     truevtx[2] = o[2];
 129     printf("True Vertex: (%f, %f, %f)\n",o[0],o[1],o[2]);
 130
 131     Double_t dNchdy = 0.;
 132     multiplicity=0.;
 133     Int_t nitstracks = 0;
 134     Int_t nitstracksFromStrange = 0;
 135     Int_t nitstracksFromStrange5or6 = 0;
 136     Int_t nitstracks1 = 0;
 137     Int_t nitstracks2 = 0;
 138     Int_t nitstracks3 = 0;
 139     Int_t nitstracks4 = 0;
 140     Int_t nitstracks5 = 0;
 141     Int_t nitstracks6 = 0;
 142     Int_t nitstracks5or6 = 0;
 143
 144     // loop on particles
 145     for(Int_t pa=0; pa<npart; pa++) {
 146       TParticle *part = (TParticle*)gAlice->GetMCApp()->Particle(pa);
 147       Int_t pdg = part->GetPdgCode();
 148       Int_t apdg = TMath::Abs(pdg);
 149       Double_t energy  = part->Energy();
 150       if(energy>6900.) continue; // reject incoming protons
 151       Double_t pz = part->Pz();
 152       Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
 153
 154       // consider charged pi,K,p,el,mu
 155       if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;
 156       // reject secondaries
 157       if(TMath::Sqrt((part->Vx()-o[0])*(part->Vx()-o[0])+(part->Vy()-o[1])*(part->Vy()-o[1]))>0.0010) continue;
 158       if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
 159     } // end loop on particles
 160     multiplicity=(Double_t)dNchdy;
 161
 162     printf(" dNch/dy = %f\n",dNchdy);
 163     //===============================================================
 164
 165
 166     esdTree->GetEvent(iev);
 167     Int_t ntracks = esd->GetNumberOfTracks();
 168     printf("Number of tracks in ESD %d   :   %d\n",iev,ntracks);
 169     if(ntracks==0) {
 170       ntVtxRes->Fill(ntracks,nitstracks5or6,nitstracksFromStrange,nitstracksFromStrange5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,nc);
 171       continue;
 172     }
 173
 174     // VERTEX
 175     AliESDVertex *vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
 176     //AliESDVertex *vertex = (AliESDVertex*)vertexer->FindPrimaryVertexOld(esd);
 177
 178     nc = (Int_t)vertex->GetNContributors();
 179     printf("nc = %d\n",nc);
 180
 181
 182     // count ITS tracks
 183     for(Int_t itrk=0; itrk<ntracks; itrk++) {
 184       AliESDtrack *esdtrack = (AliESDtrack*)esd->GetTrack(itrk);
 185       UInt_t status = esdtrack->GetStatus();
 186       // only tracks found also in ITS
 187       if(! (status&AliESDtrack::kITSin) ) continue;
 188       nitstracks++;
 189       Int_t npts = (Int_t)esdtrack->GetNcls(0);
 190       if(npts==6) {nitstracks6++;nitstracks5or6++;}
 191       if(npts==5) {nitstracks5++;nitstracks5or6++;}
 192       if(npts==4) nitstracks4++;
 193       if(npts==3) nitstracks3++;
 194       if(npts==2) nitstracks2++;
 195       if(npts==1) nitstracks1++;
 196       TParticle *part = (TParticle*)gAlice->GetMCApp()->Particle(TMath::Abs(esdtrack->GetLabel()));
 197       if(part->GetFirstMother()>=0) {
 198         Int_t mpdg = TMath::Abs(gAlice->GetMCApp()->Particle(part->GetFirstMother())->GetPdgCode());
 199         //cout<<TMath::Abs(esdtrack->GetLabel())<<"   "<<mpdg<<endl;
 200         if(mpdg==321||mpdg==311||mpdg==310||mpdg==3122||mpdg==3312||mpdg==3334) {
 201           nitstracksFromStrange++;
 202           if(npts>=5) nitstracksFromStrange5or6++;
 203         }
 204       }
 205     }
 206     printf("Number of kITSin tracks in ESD %d   :   %d\n",iev,nitstracks);
 207     printf("           6 pts: %d\n",nitstracks6);
 208     printf("           5 pts: %d\n",nitstracks5);
 209     printf("           4 pts: %d\n",nitstracks4);
 210     printf("           3 pts: %d\n",nitstracks3);
 211     printf("           2 pts: %d\n",nitstracks2);
 212     printf("           1 pts: %d\n",nitstracks1);
 213     printf("Number of kITSin tracks from s:   %d\n",nitstracksFromStrange);
 214
 215
 216     if(nc>=2) {
 217       vertex->GetXYZ(vtx);
 218       vertex->GetSigmaXYZ(errvtx);
 219       diffX = 10000.* (vtx[0] - truevtx[0]);
 220       diffY = 10000.* (vtx[1] - truevtx[1]);
 221       diffZ = 10000.* (vtx[2] - truevtx[2]);
 222       hdiffX->Fill(diffX);
 223       hdiffY->Fill(diffY);
 224       hdiffZ->Fill(diffZ);
 225       diffXerr = diffX/errvtx[0]/10000.;
 226       diffYerr = diffY/errvtx[1]/10000.;
 227       diffZerr = diffZ/errvtx[2]/10000.;
 228       hdiffXerr->Fill(diffXerr);
 229       hdiffYerr->Fill(diffYerr);
 230       hdiffZerr->Fill(diffZerr);
 231     }
 232
 233     ntVtxRes->Fill(nitstracks,nitstracks5or6,nitstracksFromStrange,nitstracksFromStrange5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,(Float_t)nc);
 234
 235   } // END LOOP ON EVENTS
 236
 237
 238   delete esdTree;
 239   delete fin;
 240   delete vertexer;
 241   delete initVertex;
 242
 243   // Output file
 244   TFile *fout  = new TFile(outname.Data(),"recreate");
 245   ntVtxRes->Write();
 246   hdiffX->Write();
 247   hdiffY->Write();
 248   hdiffZ->Write();
 249   hdiffXerr->Write();
 250   hdiffYerr->Write();
 251   hdiffZerr->Write();
 252   fout->Close();
 253
 254   return;
 255 }
 256 //----------------------------------------------------------------------------
 257 Double_t FitFunc(Double_t *x,Double_t *par) {
 258   return par[0]+par[1]/TMath::Sqrt(x[0]);
 259 }
 260 Int_t GetBin(Float_t mult) {
 261   if(mult<5.)  return 0;
 262   if(mult<7.)  return 1;
 263   if(mult<12.) return 2;
 264   if(mult<15.) return 3;
 265   if(mult<22.) return 4;
 266   return 5;
 267 }
 268 void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
 269   //---------------------------------------------------------------------
 270   // Plot efficiency, resolutions, pulls
 271   // [at least 10k pp events are needed]
 272   //---------------------------------------------------------------------
 273
 274   gStyle->SetOptStat(0);
 275   gStyle->SetOptFit(10001);
 276
 277   Int_t   nEvVtx=0;
 278   Int_t   nEvNoVtx=0;
 279   Int_t   ev,nUsedTrks;
 280   Float_t nTrks,nTrksFromStrange,nTrks5or6,nTrksFromStrange5or6,dNchdy;
 281   Float_t diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr;
 282
 283   //
 284   // HISTOGRAMS
 285   //
 286   TH1F *hdx = new TH1F("hdx","",50,-600,600);
 287   hdx->SetXTitle("#Delta x [#mu m]");
 288   hdx->SetYTitle("events");
 289   //
 290   TH1F *hdy = new TH1F("hdy","",50,-600,600);
 291   hdy->SetXTitle("#Delta y [#mu m]");
 292   hdy->SetYTitle("events");
 293   //
 294   TH1F *hdz = new TH1F("hdz","",50,-600,600);
 295   hdz->SetXTitle("#Delta z [#mu m]");
 296   hdz->SetYTitle("events");
 297   //
 298   TH1F *hTracksVtx = new TH1F("hTracksVtx","events w/ vertex",51,-0.5,50.5);
 299   hTracksVtx->SetXTitle("Number of reco tracks in TPC&ITS");
 300   hTracksVtx->SetYTitle("events");
 301   //
 302   TH1F *hTracksNoVtx = new TH1F("hTracksNoVtx","events w/o vertex",51,-0.5,50.5);
 303   hTracksNoVtx->SetXTitle("Number of reco tracks in TPC&ITS");
 304   hTracksNoVtx->SetYTitle("events");
 305   //
 306   TH1F *hTracksVtx5or6 = new TH1F("hTracksVtx5or6","events w/ vertex",51,-0.5,50.5);
 307   hTracksVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(pts>=5)");
 308   hTracksVtx5or6->SetYTitle("events");
 309   //
 310   TH1F *hTracksNoVtx5or6 = new TH1F("hTracksNoVtx5or6","events w/o vertex",51,-0.5,50.5);
 311   hTracksNoVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(pts>=5)");
 312   hTracksNoVtx5or6->SetYTitle("events");
 313   //
 314   TH1F *hTracksVtx5or6nonS = new TH1F("hTracksVtx5or6nonS","events w/ vertex",51,-0.5,50.5);
 315   hTracksVtx5or6nonS->SetXTitle("Number of reco tracks in TPC&ITS(pts>=5) NOT from s");
 316   hTracksVtx5or6nonS->SetYTitle("events");
 317   //
 318   TH1F *hTracksNoVtx5or6nonS = new TH1F("hTracksNoVtx5or6nonS","events w/o vertex",51,-0.5,50.5);
 319   hTracksNoVtx5or6nonS->SetXTitle("Number of reco tracks in TPC&ITS(pts>=5) NOT from s");
 320   hTracksNoVtx5or6nonS->SetYTitle("events");
 321   //
 322   TH1F *hPullsx = new TH1F("hPullsx","",50,-7.,7.);
 323   hPullsx->SetXTitle("#Delta x/#sqrt{<x,x>}");
 324   hPullsx->SetYTitle("events");
 325   //
 326   TH1F *hPullsy = new TH1F("hPullsy","",50,-7.,7.);
 327   hPullsy->SetXTitle("#Delta y/#sqrt{<y,y>}");
 328   hPullsy->SetYTitle("events");
 329   //
 330   TH1F *hPullsz = new TH1F("hPullsz","",50,-7.,7.);
 331   hPullsz->SetXTitle("#Delta z/#sqrt{<z,z>}");
 332   hPullsz->SetYTitle("events");
 333   //
 334
 335   Double_t mult[6]={0,0,0,0,0,0};
 336   Double_t mult2[6]={0,0,0,0,0,0};
 337   Double_t sigmamult[6]={0,0,0,0,0,0};
 338   Double_t sigmamult2[6]={0,0,0,0,0,0};
 339   Double_t tottracks[6]={0,0,0,0,0,0};
 340   Double_t fittrks[6]={0,0,0,0,0,0};
 341   Double_t tracks[6]={0,0,0,0,0,0};
 342   Double_t etracks[6]={0,0,0,0,0,0};
 343   Double_t xres[6]={0,0,0,0,0,0};
 344   Double_t exres[6]={0,0,0,0,0,0};
 345   Double_t yres[6]={0,0,0,0,0,0};
 346   Double_t eyres[6]={0,0,0,0,0,0};
 347   Double_t zres[6]={0,0,0,0,0,0};
 348   Double_t ezres[6]={0,0,0,0,0,0};
 349   Double_t xpull[6]={0,0,0,0,0,0};
 350   Double_t expull[6]={0,0,0,0,0,0};
 351   Double_t ypull[6]={0,0,0,0,0,0};
 352   Double_t eypull[6]={0,0,0,0,0,0};
 353   Double_t zpull[6]={0,0,0,0,0,0};
 354   Double_t ezpull[6]={0,0,0,0,0,0};
 355   Int_t    evts[6]={0,0,0,0,0,0};
 356   Int_t    totevts[6]={0,0,0,0,0,0};
 357   Int_t    vtx[6]={0,0,0,0,0,0};
 358   Int_t    all[6]={0,0,0,0,0,0};
 359   Double_t probVtx[6]={0,0,0,0,0,0};
 360   Double_t eprobVtx[6]={0,0,0,0,0,0};
 361   Int_t    bin;
 362   //
 363   // OTHER HISTOGRAMS
 364   //
 365   TH1F *hdx0 = new TH1F("hdx0","x resolution - bin 0",50,-500,500);
 366   TH1F *hdx1 = new TH1F("hdx1","x resolution - bin 1",50,-500,500);
 367   TH1F *hdx2 = new TH1F("hdx2","x resolution - bin 2",50,-500,500);
 368   TH1F *hdx3 = new TH1F("hdx3","x resolution - bin 3",50,-400,400);
 369   TH1F *hdx4 = new TH1F("hdx4","x resolution - bin 4",50,-300,300);
 370   TH1F *hdx5 = new TH1F("hdx5","x resolution - bin 5",50,-300,300);
 371   TH1F *hdy0 = new TH1F("hdy0","y resolution - bin 0",50,-500,500);
 372   TH1F *hdy1 = new TH1F("hdy1","y resolution - bin 1",50,-500,500);
 373   TH1F *hdy2 = new TH1F("hdy2","y resolution - bin 2",50,-500,500);
 374   TH1F *hdy3 = new TH1F("hdy3","y resolution - bin 3",50,-400,400);
 375   TH1F *hdy4 = new TH1F("hdy4","y resolution - bin 4",50,-300,300);
 376   TH1F *hdy5 = new TH1F("hdy5","y resolution - bin 5",50,-300,300);
 377   TH1F *hdz0 = new TH1F("hdz0","z resolution - bin 0",50,-1000,1000);
 378   TH1F *hdz1 = new TH1F("hdz1","z resolution - bin 1",50,-800,800);
 379   TH1F *hdz2 = new TH1F("hdz2","z resolution - bin 2",50,-800,800);
 380   TH1F *hdz3 = new TH1F("hdz3","z resolution - bin 3",50,-600,600);
 381   TH1F *hdz4 = new TH1F("hdz4","z resolution - bin 4",50,-500,500);
 382   TH1F *hdz5 = new TH1F("hdz5","z resolution - bin 5",50,-500,500);
 383
 384   TH1F *hpx0 = new TH1F("hpx0","x pull - bin 0",50,-7,7);
 385   TH1F *hpx1 = new TH1F("hpx1","x pull - bin 1",50,-7,7);
 386   TH1F *hpx2 = new TH1F("hpx2","x pull - bin 2",50,-7,7);
 387   TH1F *hpx3 = new TH1F("hpx3","x pull - bin 3",50,-7,7);
 388   TH1F *hpx4 = new TH1F("hpx4","x pull - bin 4",50,-7,7);
 389   TH1F *hpx5 = new TH1F("hpx5","x pull - bin 5",50,-7,7);
 390   TH1F *hpy0 = new TH1F("hpy0","y pull - bin 0",50,-7,7);
 391   TH1F *hpy1 = new TH1F("hpy1","y pull - bin 1",50,-7,7);
 392   TH1F *hpy2 = new TH1F("hpy2","y pull - bin 2",50,-7,7);
 393   TH1F *hpy3 = new TH1F("hpy3","y pull - bin 3",50,-7,7);
 394   TH1F *hpy4 = new TH1F("hpy4","y pull - bin 4",50,-7,7);
 395   TH1F *hpy5 = new TH1F("hpy5","y pull - bin 5",50,-7,7);
 396   TH1F *hpz0 = new TH1F("hpz0","z pull - bin 0",50,-7,7);
 397   TH1F *hpz1 = new TH1F("hpz1","z pull - bin 1",50,-7,7);
 398   TH1F *hpz2 = new TH1F("hpz2","z pull - bin 2",50,-7,7);
 399   TH1F *hpz3 = new TH1F("hpz3","z pull - bin 3",50,-7,7);
 400   TH1F *hpz4 = new TH1F("hpz4","z pull - bin 4",50,-7,7);
 401   TH1F *hpz5 = new TH1F("hpz5","z pull - bin 5",50,-7,7);
 402
 403
 404   TH1F *hmult = new TH1F("hmult","hmult",50,-0.5,49.5);
 405
 406   TFile *in = new TFile(inname);
 407   TNtuple *nt = (TNtuple*)in->Get("ntVtxRes");
 408   nt->SetBranchAddress("ntracks",&nTrks);
 409   nt->SetBranchAddress("nitstracks5or6",&nTrks5or6);
 410   nt->SetBranchAddress("nitstracksFromStrange",&nTrksFromStrange);
 411   nt->SetBranchAddress("nitstracksFromStrange5or6",&nTrksFromStrange5or6);
 412   nt->SetBranchAddress("diffX",&diffX);
 413   nt->SetBranchAddress("diffY",&diffY);
 414   nt->SetBranchAddress("diffZ",&diffZ);
 415   nt->SetBranchAddress("diffXerr",&diffXerr);
 416   nt->SetBranchAddress("diffYerr",&diffYerr);
 417   nt->SetBranchAddress("diffZerr",&diffZerr);
 418   nt->SetBranchAddress("multiplicity",&dNchdy);
 419   Int_t entries = (Int_t)nt->GetEntries();
 420   Int_t nbytes=0;
 421
 422   for(Int_t i=0; i<entries; i++) {
 423     nbytes += nt->GetEvent(i);
 424     bin = GetBin(dNchdy);
 425     mult[bin] += dNchdy;
 426     hmult->Fill(dNchdy);
 427     totevts[bin]++;
 428     if(diffX < 90000.) { // vtx found
 429       nEvVtx++;
 430       mult2[bin] += dNchdy;
 431       vtx[bin]++;
 432       tottracks[bin] += nTrks;
 433       evts[bin]++;
 434       if(bin==0) hdx0->Fill(diffX);
 435       if(bin==1) hdx1->Fill(diffX);
 436       if(bin==2) hdx2->Fill(diffX);
 437       if(bin==3) hdx3->Fill(diffX);
 438       if(bin==4) hdx4->Fill(diffX);
 439       if(bin==5) hdx5->Fill(diffX);
 440       if(bin==0) hdy0->Fill(diffY);
 441       if(bin==1) hdy1->Fill(diffY);
 442       if(bin==2) hdy2->Fill(diffY);
 443       if(bin==3) hdy3->Fill(diffY);
 444       if(bin==4) hdy4->Fill(diffY);
 445       if(bin==5) hdy5->Fill(diffY);
 446       if(bin==0) hdz0->Fill(diffZ);
 447       if(bin==1) hdz1->Fill(diffZ);
 448       if(bin==2) hdz2->Fill(diffZ);
 449       if(bin==3) hdz3->Fill(diffZ);
 450       if(bin==4) hdz4->Fill(diffZ);
 451       if(bin==5) hdz5->Fill(diffZ);
 452       hdx->Fill(diffX);
 453       hdy->Fill(diffY);
 454       hdz->Fill(diffZ);
 455       hTracksVtx->Fill(nTrks);
 456       hTracksVtx5or6->Fill(nTrks5or6);
 457       hTracksVtx5or6nonS->Fill(nTrks5or6-nTrksFromStrange5or6);
 458       hPullsx->Fill(diffXerr);
 459       hPullsy->Fill(diffYerr);
 460       hPullsz->Fill(diffZerr);
 461       if(bin==0) hpx0->Fill(diffXerr);
 462       if(bin==1) hpx1->Fill(diffXerr);
 463       if(bin==2) hpx2->Fill(diffXerr);
 464       if(bin==3) hpx3->Fill(diffXerr);
 465       if(bin==4) hpx4->Fill(diffXerr);
 466       if(bin==5) hpx5->Fill(diffXerr);
 467       if(bin==0) hpy0->Fill(diffYerr);
 468       if(bin==1) hpy1->Fill(diffYerr);
 469       if(bin==2) hpy2->Fill(diffYerr);
 470       if(bin==3) hpy3->Fill(diffYerr);
 471       if(bin==4) hpy4->Fill(diffYerr);
 472       if(bin==5) hpy5->Fill(diffYerr);
 473       if(bin==0) hpz0->Fill(diffZerr);
 474       if(bin==1) hpz1->Fill(diffZerr);
 475       if(bin==2) hpz2->Fill(diffZerr);
 476       if(bin==3) hpz3->Fill(diffZerr);
 477       if(bin==4) hpz4->Fill(diffZerr);
 478       if(bin==5) hpz5->Fill(diffZerr);
 479     } else {
 480       nEvNoVtx++;
 481       hTracksNoVtx->Fill(nTrks);
 482       hTracksNoVtx5or6->Fill(nTrks5or6);
 483       hTracksNoVtx5or6nonS->Fill(nTrks5or6-nTrksFromStrange5or6);
 484     }
 485     all[bin]++;
 486   }
 487
 488
 489   for(Int_t k=0;k<6;k++) {
 490     mult2[k]      /= evts[k];
 491     mult[k]      /= totevts[k];
 492     tottracks[k] /= evts[k];
 493     probVtx[k]    = (Double_t)vtx[k]/all[k];
 494     eprobVtx[k]   = (Double_t)vtx[k]/all[k]/all[k]+(Double_t)vtx[k]*vtx[k]/all[k]/all[k]/all[k];
 495     eprobVtx[k] = TMath::Sqrt(eprobVtx[k]);
 496   }
 497
 498   // calculate spread in multiplicity
 499   for(Int_t i=0; i<entries; i++) {
 500     nbytes += nt->GetEvent(i);
 501     bin = GetBin(dNchdy);
 502     sigmamult[bin] += (dNchdy-mult[bin])*(dNchdy-mult[bin])/totevts[bin];
 503     if(diffX < 90000.) sigmamult2[bin] += (dNchdy-mult2[bin])*(dNchdy-mult2[bin])/evts[bin];
 504   }
 505
 506   for(Int_t k=0;k<6;k++) {
 507     sigmamult[k]  = TMath::Sqrt(sigmamult[k]);
 508     sigmamult2[k] = TMath::Sqrt(sigmamult2[k]);
 509   }
 510
 511   // fraction of found vertices
 512   Float_t frac = (Float_t)nEvVtx/(nEvVtx+nEvNoVtx);
 513   printf(" Fraction of events with vertex: %f\n",frac);
 514
 515
 516   // FIT FUNCTIONS
 517   TF1 *g = new TF1("g","gaus");
 518   g->SetLineWidth(3);
 519   TF1 *line = new TF1("line","pol1");
 520   line->SetLineWidth(3);
 521   TF1 *func = new TF1("func",FitFunc,2,35,2);
 522   func->SetLineWidth(1);
 523
 524   Double_t x1,y1,sigmafit;
 525   Char_t message[100];
 526
 527   //
 528   //            DRAW PLOTS
 529   //
 530   gStyle->SetOptFit(1111);
 531
 532   // Tracks in ITS for events w/ and w/o vertex
 533   TCanvas *c1 = new TCanvas("c1","c1",0,0,1000,500);
 534   c1->Divide(2,1);
 535   c1->cd(1);
 536   hTracksVtx->Draw();
 537   c1->cd(2);
 538   hTracksNoVtx->Draw();
 539
 540   // probability vs ntracks
 541   TCanvas *c1b = new TCanvas("c1b","c1b",0,0,500,500);
 542   TH1F *hTotTracks = (TH1F*)hTracksVtx->Clone("hTotTracks");
 543   hTotTracks->Add(hTracksNoVtx);
 544   TH1F *hProbTracks = (TH1F*)hTracksVtx->Clone("hProbTracks");
 545   hProbTracks->Divide(hTotTracks);
 546   hProbTracks->SetYTitle("Probability");
 547   hProbTracks->SetTitle("Probability to find the vertex");
 548   hProbTracks->Draw();
 549
 550   TCanvas *c1c = new TCanvas("c1c","c1c",0,0,500,500);
 551   TH1F *hTotTracks5or6 = (TH1F*)hTracksVtx5or6->Clone("hTotTracks5or6");
 552   hTotTracks5or6->Add(hTracksNoVtx5or6);
 553   TH1F *hProbTracks5or6 = (TH1F*)hTracksVtx5or6->Clone("hProbTracks5or6");
 554   hProbTracks5or6->Divide(hTotTracks5or6);
 555   hProbTracks5or6->SetYTitle("Probability");
 556   hProbTracks5or6->SetTitle("Probability to find the vertex");
 557   hProbTracks5or6->Draw();
 558
 559   TCanvas *c1d = new TCanvas("c1d","c1d",0,0,500,500);
 560   TH1F *hTotTracks5or6nonS = (TH1F*)hTracksVtx5or6nonS->Clone("hTotTracks5or6nonS");
 561   hTotTracks5or6nonS->Add(hTracksNoVtx5or6nonS);
 562   TH1F *hProbTracks5or6nonS = (TH1F*)hTracksVtx5or6nonS->Clone("hProbTracks5or6nonS");
 563   hProbTracks5or6nonS->Divide(hTotTracks5or6nonS);
 564   hProbTracks5or6nonS->SetYTitle("Probability");
 565   hProbTracks5or6nonS->SetTitle("Probability to find the vertex");
 566   hProbTracks5or6nonS->Draw();
 567
 568   // Global resolutions
 569   TCanvas *c2 = new TCanvas("c2","c2",0,0,1000,400);
 570   c2->Divide(3,1);
 571   c2->cd(1);
 572   hdx->Draw();
 573   g->SetRange(-1.*hdx->GetRMS(),+1.*hdx->GetRMS());
 574   hdx->Fit("g","R");
 575   c2->cd(2);
 576   hdy->Draw();
 577   g->SetRange(-1.*hdy->GetRMS(),+1.*hdy->GetRMS());
 578   hdy->Fit("g","R");
 579   c2->cd(3);
 580   hdz->Draw();
 581   g->SetRange(-1.*hdz->GetRMS(),+1.*hdz->GetRMS());
 582   hdz->Fit("g","R");
 583
 584   // Pulls
 585   TCanvas *c4 = new TCanvas("c4","c4",0,0,1000,400);
 586   c4->Divide(3,1);
 587   c4->cd(1);
 588   hPullsx->Draw();
 589   g->SetRange(-3.*hPullsx->GetRMS(),+3.*hPullsx->GetRMS());
 590   hPullsx->Fit("g","R");
 591   c4->cd(2);
 592   hPullsy->Draw();
 593   g->SetRange(-3.*hPullsy->GetRMS(),+3.*hPullsy->GetRMS());
 594   hPullsy->Fit("g","R");
 595   c4->cd(3);
 596   hPullsz->Draw();
 597   g->SetRange(-3.*hPullsz->GetRMS(),+3.*hPullsz->GetRMS());
 598   hPullsz->Fit("g","R");
 599
 600   // probability VS multiplicity
 601   TCanvas *c5 = new TCanvas("c5","c5",0,0,600,500);
 602   TH2F *fr5 = new TH2F("fr5","",2,0,35,2,0,1.1);
 603   fr5->SetXTitle("dN_{ch}/dy");
 604   fr5->SetYTitle("Probability");
 605   fr5->Draw();
 606   TGraphErrors *gr5 = new TGraphErrors(6,mult,probVtx,sigmamult,eprobVtx);
 607   gr5->Draw("p");
 608   gr5->SetMarkerStyle(21);
 609
 610   //
 611   // resolutions VS multiplicity
 612   //
 613   // X
 614   TCanvas *ax = new TCanvas("ax","ax",0,0,900,700);
 615   ax->Divide(3,2);
 616   ax->cd(1);
 617   hdx0->Draw();
 618   g->SetRange(-3.*hdx0->GetRMS(),+3.*hdx0->GetRMS());
 619   hdx0->Fit("g","Q");
 620   xres[0]=g->GetParameter(2);
 621   exres[0]=g->GetParError(2);
 622   ax->cd(2);
 623   hdx1->Draw();
 624   g->SetRange(-3.*hdx1->GetRMS(),+3.*hdx1->GetRMS());
 625   hdx1->Fit("g","Q");
 626   xres[1]=g->GetParameter(2);
 627   exres[1]=g->GetParError(2);
 628   ax->cd(3);
 629   hdx2->Draw();
 630   g->SetRange(-3.*hdx2->GetRMS(),+3.*hdx2->GetRMS());
 631   hdx2->Fit("g","Q");
 632   xres[2]=g->GetParameter(2);
 633   exres[2]=g->GetParError(2);
 634   ax->cd(4);
 635   hdx3->Draw();
 636   g->SetRange(-3.*hdx3->GetRMS(),+3.*hdx3->GetRMS());
 637   hdx3->Fit("g","Q");
 638   xres[3]=g->GetParameter(2);
 639   exres[3]=g->GetParError(2);
 640   ax->cd(5);
 641   hdx4->Draw();
 642   g->SetRange(-3.*hdx4->GetRMS(),+3.*hdx4->GetRMS());
 643   hdx4->Fit("g","Q");
 644   xres[4]=g->GetParameter(2);
 645   exres[4]=g->GetParError(2);
 646   ax->cd(6);
 647   hdx5->Draw();
 648   g->SetRange(-3.*hdx5->GetRMS(),+3.*hdx5->GetRMS());
 649   hdx5->Fit("g","Q");
 650   xres[5]=g->GetParameter(2);
 651   exres[5]=g->GetParError(2);
 652   // Y
 653   TCanvas *ay = new TCanvas("ay","ay",0,0,900,700);
 654   ay->Divide(3,2);
 655   ay->cd(1);
 656   hdy0->Draw();
 657   g->SetRange(-3.*hdy0->GetRMS(),+3.*hdy0->GetRMS());
 658   hdy0->Fit("g","Q");
 659   yres[0]=g->GetParameter(2);
 660   eyres[0]=g->GetParError(2);
 661   ay->cd(2);
 662   hdy1->Draw();
 663   g->SetRange(-3.*hdy1->GetRMS(),+3.*hdy1->GetRMS());
 664   hdy1->Fit("g","Q");
 665   yres[1]=g->GetParameter(2);
 666   eyres[1]=g->GetParError(2);
 667   ay->cd(3);
 668   hdy2->Draw();
 669   g->SetRange(-3.*hdy2->GetRMS(),+3.*hdy2->GetRMS());
 670   hdy2->Fit("g","Q");
 671   yres[2]=g->GetParameter(2);
 672   eyres[2]=g->GetParError(2);
 673   ay->cd(4);
 674   hdy3->Draw();
 675   g->SetRange(-3.*hdy3->GetRMS(),+3.*hdy3->GetRMS());
 676   hdy3->Fit("g","Q");
 677   yres[3]=g->GetParameter(2);
 678   eyres[3]=g->GetParError(2);
 679   ay->cd(5);
 680   hdy4->Draw();
 681   g->SetRange(-3.*hdy4->GetRMS(),+3.*hdy4->GetRMS());
 682   hdy4->Fit("g","Q");
 683   yres[4]=g->GetParameter(2);
 684   eyres[4]=g->GetParError(2);
 685   ay->cd(6);
 686   hdy5->Draw();
 687   g->SetRange(-3.*hdy5->GetRMS(),+3.*hdy5->GetRMS());
 688   hdy5->Fit("g","Q");
 689   yres[5]=g->GetParameter(2);
 690   eyres[5]=g->GetParError(2);
 691   // Z
 692   TCanvas *az = new TCanvas("az","az",0,0,900,700);
 693   az->Divide(3,2);
 694   az->cd(1);
 695   hdz0->Draw();
 696   g->SetRange(-3.*hdz0->GetRMS(),+3.*hdz0->GetRMS());
 697   hdz0->Fit("g","Q");
 698   zres[0]=g->GetParameter(2);
 699   ezres[0]=g->GetParError(2);
 700   az->cd(2);
 701   hdz1->Draw();
 702   g->SetRange(-3.*hdz1->GetRMS(),+3.*hdz1->GetRMS());
 703   hdz1->Fit("g","Q");
 704   zres[1]=g->GetParameter(2);
 705   ezres[1]=g->GetParError(2);
 706   az->cd(3);
 707   hdz2->Draw();
 708   g->SetRange(-3.*hdz2->GetRMS(),+3.*hdz2->GetRMS());
 709   hdz2->Fit("g","Q");
 710   zres[2]=g->GetParameter(2);
 711   ezres[2]=g->GetParError(2);
 712   az->cd(4);
 713   hdz3->Draw();
 714   g->SetRange(-3.*hdz3->GetRMS(),+3.*hdz3->GetRMS());
 715   hdz3->Fit("g","Q");
 716   zres[3]=g->GetParameter(2);
 717   ezres[3]=g->GetParError(2);
 718   az->cd(5);
 719   hdz4->Draw();
 720   g->SetRange(-3.*hdz4->GetRMS(),+3.*hdz4->GetRMS());
 721   hdz4->Fit("g","Q");
 722   zres[4]=g->GetParameter(2);
 723   ezres[4]=g->GetParError(2);
 724   az->cd(6);
 725   hdz5->Draw();
 726   g->SetRange(-3.*hdz5->GetRMS(),+3.*hdz5->GetRMS());
 727   hdz5->Fit("g","Q");
 728   zres[5]=g->GetParameter(2);
 729   ezres[5]=g->GetParError(2);
 730
 731
 732   //
 733   // pulls VS multiplicity
 734   //
 735   // X
 736   TCanvas *bx = new TCanvas("bx","bx",0,0,900,700);
 737   bx->Divide(3,2);
 738   bx->cd(1);
 739   hpx0->Draw();
 740   g->SetRange(-3.*hpx0->GetRMS(),+3.*hpx0->GetRMS());
 741   hpx0->Fit("g","Q");
 742   xpull[0]=g->GetParameter(2);
 743   expull[0]=g->GetParError(2);
 744   bx->cd(2);
 745   hpx1->Draw();
 746   g->SetRange(-3.*hpx1->GetRMS(),+3.*hpx1->GetRMS());
 747   hpx1->Fit("g","Q");
 748   xpull[1]=g->GetParameter(2);
 749   expull[1]=g->GetParError(2);
 750   bx->cd(3);
 751   hpx2->Draw();
 752   g->SetRange(-3.*hpx2->GetRMS(),+3.*hpx2->GetRMS());
 753   hpx2->Fit("g","Q");
 754   xpull[2]=g->GetParameter(2);
 755   expull[2]=g->GetParError(2);
 756   bx->cd(4);
 757   hpx3->Draw();
 758   g->SetRange(-3.*hpx3->GetRMS(),+3.*hpx3->GetRMS());
 759   hpx3->Fit("g","Q");
 760   xpull[3]=g->GetParameter(2);
 761   expull[3]=g->GetParError(2);
 762   bx->cd(5);
 763   hpx4->Draw();
 764   g->SetRange(-3.*hpx4->GetRMS(),+3.*hpx4->GetRMS());
 765   hpx4->Fit("g","Q");
 766   xpull[4]=g->GetParameter(2);
 767   expull[4]=g->GetParError(2);
 768   bx->cd(6);
 769   hpx5->Draw();
 770   g->SetRange(-3.*hpx5->GetRMS(),+3.*hpx5->GetRMS());
 771   hpx5->Fit("g","Q");
 772   xpull[5]=g->GetParameter(2);
 773   expull[5]=g->GetParError(2);
 774   // Y
 775   TCanvas *by = new TCanvas("by","by",0,0,900,700);
 776   by->Divide(3,2);
 777   by->cd(1);
 778   hpy0->Draw();
 779   g->SetRange(-3.*hpy0->GetRMS(),+3.*hpy0->GetRMS());
 780   hpy0->Fit("g","Q");
 781   ypull[0]=g->GetParameter(2);
 782   eypull[0]=g->GetParError(2);
 783   by->cd(2);
 784   hpy1->Draw();
 785   g->SetRange(-3.*hpy1->GetRMS(),+3.*hpy1->GetRMS());
 786   hpy1->Fit("g","Q");
 787   ypull[1]=g->GetParameter(2);
 788   eypull[1]=g->GetParError(2);
 789   by->cd(3);
 790   hpy2->Draw();
 791   g->SetRange(-3.*hpy2->GetRMS(),+3.*hpy2->GetRMS());
 792   hpy2->Fit("g","Q");
 793   ypull[2]=g->GetParameter(2);
 794   eypull[2]=g->GetParError(2);
 795   by->cd(4);
 796   hpy3->Draw();
 797   g->SetRange(-3.*hpy3->GetRMS(),+3.*hpy3->GetRMS());
 798   hpy3->Fit("g","Q");
 799   ypull[3]=g->GetParameter(2);
 800   eypull[3]=g->GetParError(2);
 801   by->cd(5);
 802   hpy4->Draw();
 803   g->SetRange(-3.*hpy4->GetRMS(),+3.*hpy4->GetRMS());
 804   hpy4->Fit("g","Q");
 805   ypull[4]=g->GetParameter(2);
 806   eypull[4]=g->GetParError(2);
 807   by->cd(6);
 808   hpy5->Draw();
 809   g->SetRange(-3.*hpy5->GetRMS(),+3.*hpy5->GetRMS());
 810   hpy5->Fit("g","Q");
 811   ypull[5]=g->GetParameter(2);
 812   eypull[5]=g->GetParError(2);
 813   // Z
 814   TCanvas *bz = new TCanvas("bz","bz",0,0,900,700);
 815   bz->Divide(3,2);
 816   bz->cd(1);
 817   hpz0->Draw();
 818   g->SetRange(-3.*hpz0->GetRMS(),+3.*hpz0->GetRMS());
 819   hpz0->Fit("g","Q");
 820   zpull[0]=g->GetParameter(2);
 821   ezpull[0]=g->GetParError(2);
 822   bz->cd(2);
 823   hpz1->Draw();
 824   g->SetRange(-3.*hpz1->GetRMS(),+3.*hpz1->GetRMS());
 825   hpz1->Fit("g","Q");
 826   zpull[1]=g->GetParameter(2);
 827   ezpull[1]=g->GetParError(2);
 828   bz->cd(3);
 829   hpz2->Draw();
 830   g->SetRange(-3.*hpz2->GetRMS(),+3.*hpz2->GetRMS());
 831   hpz2->Fit("g","Q");
 832   zpull[2]=g->GetParameter(2);
 833   ezpull[2]=g->GetParError(2);
 834   bz->cd(4);
 835   hpz3->Draw();
 836   g->SetRange(-3.*hpz3->GetRMS(),+3.*hpz3->GetRMS());
 837   hpz3->Fit("g","Q");
 838   zpull[3]=g->GetParameter(2);
 839   ezpull[3]=g->GetParError(2);
 840   bz->cd(5);
 841   hpz4->Draw();
 842   g->SetRange(-3.*hpz4->GetRMS(),+3.*hpz4->GetRMS());
 843   hpz4->Fit("g","Q");
 844   zpull[4]=g->GetParameter(2);
 845   ezpull[4]=g->GetParError(2);
 846   bz->cd(6);
 847   hpz5->Draw();
 848   g->SetRange(-3.*hpz5->GetRMS(),+3.*hpz5->GetRMS());
 849   hpz5->Fit("g","Q");
 850   zpull[5]=g->GetParameter(2);
 851   ezpull[5]=g->GetParError(2);
 852
 853   gStyle->SetOptFit(0);
 854
 855   // res VS dNchdy
 856   TCanvas *c6 = new TCanvas("c6","c6",0,0,600,500);
 857   TH2F *fr6 = new TH2F("fr6","",2,0,35,2,0,200);
 858   fr6->SetXTitle("dN_{ch}/dy");
 859   fr6->SetYTitle("#sigma [#mu m]");
 860   fr6->Draw();
 861   sigmamult2[0]=sigmamult2[1]=sigmamult2[2]=sigmamult2[3]=sigmamult2[4]=sigmamult2[5]=0.;
 862   TGraphErrors *gr6x = new TGraphErrors(6,mult2,xres,sigmamult2,exres);
 863   gr6x->Draw("p");
 864   gr6x->SetMarkerStyle(22);
 865   gr6x->Fit("func","E,R");
 866   TGraphErrors *gr6z = new TGraphErrors(6,mult2,zres,sigmamult2,ezres);
 867   gr6z->Draw("p");
 868   gr6z->SetMarkerStyle(26);
 869   gr6z->Fit("func","E,R");
 870   TLegend *leg6 = new TLegend(0.6,0.75,0.9,0.9);
 871   leg6->AddEntry(gr6x,"x/y coordinate","p");
 872   leg6->AddEntry(gr6z,"z coordinate","p");
 873   leg6->SetFillStyle(0);
 874   leg6->SetBorderSize(0);
 875   leg6->Draw();
 876
 877   // pull VS dNchdy
 878   TCanvas *c8 = new TCanvas("c8","c8",0,0,600,500);
 879   TH2F *fr8 = new TH2F("fr8","",2,0,35,2,0,2.);
 880   fr8->SetXTitle("dN_{ch}/dy");
 881   fr8->SetYTitle("pull");
 882   fr8->Draw();
 883   TGraphErrors *gr8x = new TGraphErrors(6,mult2,xpull,sigmamult2,expull);
 884   gr8x->Draw("p");
 885   gr8x->SetMarkerStyle(22);
 886   TGraphErrors *gr8z = new TGraphErrors(6,mult2,zpull,sigmamult2,ezpull);
 887   gr8z->Draw("p");
 888   gr8z->SetMarkerStyle(26);
 889   TLegend *leg8 = new TLegend(0.6,0.75,0.9,0.9);
 890   leg8->AddEntry(gr6x,"x/y coordinate","p");
 891   leg8->AddEntry(gr6z,"z coordinate","p");
 892   leg8->SetFillStyle(0);
 893   leg8->SetBorderSize(0);
 894   leg8->Draw();
 895   TLine *l8 = new TLine(0,1,35,1);
 896   l8->SetLineStyle(2);
 897   l8->Draw();
 898
 899   return;
 900 }
 901 //----------------------------------------------------------------------------
 902 void AliComputeVtxMeanFromESD(TString file="AliESDs.root",
 903                               Int_t mincontr=20) {
 904   //-----------------------------------------------------------------------
 905   // Compute weighted mean and cov. matrix from a set of AliESD
 906   //-----------------------------------------------------------------------
 907
 908   gStyle->SetOptStat(0);
 909   //gStyle->SetOptFit(0);
 910
 911
 912   Double_t vtx[3],covvtx[6];
 913   TTree *esdTree = 0;
 914   AliESD *esd = 0;
 915   AliESDVertex *vertex = 0;
 916   TString vtitle;
 917   Int_t nc,events,total=0;
 918
 919   Double_t avx=0.;
 920   Double_t avy=0.;
 921   Double_t avz=0.;
 922   Double_t wgtavx=0.;
 923   Double_t wgtavy=0.;
 924   Double_t wgtavz=0.;
 925   Double_t sum=0.;
 926   Double_t sumwgtx=0.;
 927   Double_t sumwgty=0.;
 928   Double_t sumwgtz=0.;
 929   Double_t rmsx=0;
 930   Double_t rmsy=0;
 931   Double_t rmsz=0;
 932   Double_t varx=0;
 933   Double_t vary=0;
 934   Double_t varz=0;
 935   Double_t covxy=0;
 936   Double_t covxz=0;
 937   Double_t covyz=0;
 938   Double_t eavx,eavy,eavz,ewgtavx,ewgtavy,ewgtavz;
 939
 940   TH1F* hx = new TH1F("hx","",200,-1000,1000);
 941   hx->SetXTitle("vertex x [#mu m]");
 942   hx->SetYTitle("events");
 943   TH1F* hy = new TH1F("hy","",200,-1000,1000);
 944   hy->SetXTitle("vertex y [#mu m]");
 945   hy->SetYTitle("events");
 946   TH1F* hz = new TH1F("hz","",200,-10,10);
 947   hz->SetXTitle("vertex z [cm]");
 948   hz->SetYTitle("events");
 949
 950
 951   for(Int_t i=1; i<100; i++) {
 952     TString inname="event.";
 953     inname+=i;
 954     inname.Append("/");
 955     inname.Append(file.Data());
 956     if(gSystem->AccessPathName(inname.Data())) continue;
 957     TFile *fin = TFile::Open(inname.Data());
 958     esdTree = (TTree*)fin->Get("esdTree");
 959     if(!esdTree) continue;
 960     esdTree->SetBranchAddress("ESD",&esd);
 961     events = esdTree->GetEntries();
 962     printf("Number of events in ESD tree: %d\n",events);
 963     total += events;
 964
 965     for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
 966       esdTree->GetEvent(iev);
 967       vertex = (AliESDVertex*)esd->GetPrimaryVertex();
 968       if(!vertex) continue;
 969       nc = vertex->GetNContributors();
 970       if(nc<mincontr) continue;
 971       vertex->GetXYZ(vtx);
 972       vertex->GetSigmaXYZ(errvtx);
 973       //printf("%f %f %f %f %f %f\n",vtx[0],vtx[1],vtx[2],errvtx[0],errvtx[1],errvtx[2]);
 974
 975       avx += vtx[0];
 976       avy += vtx[1];
 977       avz += vtx[2];
 978       sum += 1.;
 979       wgtavx += vtx[0]/errvtx[0]/errvtx[0];
 980       wgtavy += vtx[1]/errvtx[1]/errvtx[1];
 981       wgtavz += vtx[2]/errvtx[2]/errvtx[2];
 982       sumwgtx += 1./errvtx[0]/errvtx[0];
 983       sumwgty += 1./errvtx[1]/errvtx[1];
 984       sumwgtz += 1./errvtx[2]/errvtx[2];
 985
 986       hx->Fill(1e4*vtx[0]);
 987       hy->Fill(1e4*vtx[1]);
 988       hz->Fill(vtx[2]);
 989
 990       vertex = 0;
 991     }
 992     esdTree=0;
 993     fin->Close();
 994     delete fin;
 995   }
 996
 997   avx /= sum;
 998   avy /= sum;
 999   avz /= sum;
1000   wgtavx /= sumwgtx;
1001   wgtavy /= sumwgty;
1002   wgtavz /= sumwgtz;
1003   ewgtavx = 1./TMath::Sqrt(sumwgtx);
1004   ewgtavy = 1./TMath::Sqrt(sumwgty);
1005   ewgtavz = 1./TMath::Sqrt(sumwgtz);
1006
1007
1008   for(Int_t i=1; i<100; i++) {
1009     TString inname="event.";
1010     inname+=i;
1011     inname.Append("/");
1012     inname.Append(file.Data());
1013     if(gSystem->AccessPathName(inname.Data())) continue;
1014     TFile *fin = TFile::Open(inname.Data());
1015     esdTree = (TTree*)fin->Get("esdTree");
1016     if(!esdTree) continue;
1017     esdTree->SetBranchAddress("ESD",&esd);
1018     events = esdTree->GetEntries();
1019     for(Int_t iev=0; iev<events; iev++) { //LOOP ON EVENTS
1020       esdTree->GetEvent(iev);
1021       vertex = (AliESDVertex*)esd->GetPrimaryVertex();
1022       if(!vertex) continue;
1023       nc = vertex->GetNContributors();
1024       if(nc<mincontr) continue;
1025       vertex->GetXYZ(vtx);
1026       vertex->GetSigmaXYZ(errvtx);
1027
1028       varx += (vtx[0]-avx)*(vtx[0]-avx);
1029       vary += (vtx[1]-avy)*(vtx[1]-avy);
1030       varz += (vtx[2]-avz)*(vtx[2]-avz);
1031       covxy += (vtx[0]-avx)*(vtx[1]-avy);
1032       covxz += (vtx[0]-avx)*(vtx[2]-avz);
1033       covyz += (vtx[1]-avy)*(vtx[2]-avz);
1034
1035       vertex = 0;
1036     }
1037   }
1038
1039   varx /= sum;
1040   vary /= sum;
1041   varz /= sum;
1042   covxy /= sum;
1043   covxz /= sum;
1044   covyz /= sum;
1045   rmsx = TMath::Sqrt(varx);
1046   rmsy = TMath::Sqrt(vary);
1047   rmsz = TMath::Sqrt(varz);
1048   eavx = rmsx/TMath::Sqrt(sum);
1049   eavy = rmsy/TMath::Sqrt(sum);
1050   eavz = rmsz/TMath::Sqrt(sum);
1051
1052
1053   printf("\n\nNumber of events: Total %d, Used %d\n",total,sum);
1054   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);
1055   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);
1056   printf("RMS:\n x = %f mum\n y = %f mum\n z = %f mum\n",1e4*rmsx,1e4*rmsy,1e4*rmsz);
1057
1058   TCanvas *c = new TCanvas("c","c",0,0,1000,500);
1059   c->Divide(3,1);
1060   c->cd(1);
1061   hx->Draw();
1062   TF1 *gx = new TF1("gx","gaus",-1000,1000);
1063   gx->SetLineColor(2);
1064   hx->Fit(gx);
1065   TF1 *gxx = (TF1*)gx->Clone("gxx");
1066   gxx->FixParameter(2,50.);
1067   gxx->SetLineStyle(2);
1068   gxx->Draw("same");
1069   c->cd(2);
1070   hy->Draw();
1071   TF1 *gy = new TF1("gy","gaus",-1000,1000);
1072   gy->SetLineColor(2);
1073   hy->Fit(gy);
1074   TF1 *gyy = (TF1*)gy->Clone("gyy");
1075   gyy->FixParameter(2,50.);
1076   gyy->SetLineStyle(2);
1077   gyy->Draw("same");
1078   c->cd(3);
1079   hz->Draw();
1080   TF1 *gz = new TF1("gz","gaus",-10,10);
1081   gz->SetLineColor(2);
1082   hz->Fit(gz);
1083   TF1 *gzz = (TF1*)gz->Clone("gzz");
1084   gzz->FixParameter(2,5.3);
1085   gzz->SetLineStyle(2);
1086   gzz->Draw("same");
1087
1088   vtx[0] = wgtavx;
1089   vtx[1] = wgtavy;
1090   vtx[2] = wgtavz;
1091   covvtx[0] = varx;
1092   covvtx[1] = covxy;
1093   covvtx[2] = vary;
1094   covvtx[3] = covxz;
1095   covvtx[4] = covyz;
1096   covvtx[5] = 5.3;
1097   AliESDVertex *vtxmean = new AliESDVertex(vtx,covvtx,0.,(Int_t)sum);
1098   vtxmean->SetTitle("Vertex from weighted average");
1099   TFile *out = new TFile("AliESDVertexMean.root","recreate");
1100   vtxmean->Write("vtxmean");
1101   out->Close();
1102
1103   return;
1104 }
1105 //----------------------------------------------------------------------------
1106 void AliComputeVtxMeanFromTree(TString file="AliESDVertices.root",
1107                                Int_t mincontr=20) {
1108   //-----------------------------------------------------------------------
1109   // Compute weighted mean and cov. matrix from a set of AliESDVertices
1110   //-----------------------------------------------------------------------
1111   gStyle->SetOptStat(0);
1112   //gStyle->SetOptFit(0);
1113
1114
1115   Double_t vtx[3],covvtx[6];
1116   TString vtitle;
1117   Int_t nc,total=0;
1118
1119   Double_t avx=0.;
1120   Double_t avy=0.;
1121   Double_t avz=0.;
1122   Double_t wgtavx=0.;
1123   Double_t wgtavy=0.;
1124   Double_t wgtavz=0.;
1125   Double_t sum=0.;
1126   Double_t sumwgtx=0.;
1127   Double_t sumwgty=0.;
1128   Double_t sumwgtz=0.;
1129   Double_t rmsx=0;
1130   Double_t rmsy=0;
1131   Double_t rmsz=0;
1132   Double_t eavx,eavy,eavz,ewgtavx,ewgtavy,ewgtavz;
1133
1134   TH1F* hx = new TH1F("hx","",200,-1000,1000);
1135   hx->SetXTitle("vertex x [#mu m]");
1136   hx->SetYTitle("events");
1137   TH1F* hy = new TH1F("hy","",200,-1000,1000);
1138   hy->SetXTitle("vertex y [#mu m]");
1139   hy->SetYTitle("events");
1140   TH1F* hz = new TH1F("hz","",200,-10,10);
1141   hz->SetXTitle("vertex z [cm]");
1142   hz->SetYTitle("events");
1143
1144
1145   if(gSystem->AccessPathName(file.Data())) continue;
1146   TFile *fin = TFile::Open(file.Data());
1147   TTree *tree = (TTree*)fin->Get("TreeV");
1148   AliESDVertex *vertex = 0;
1149   tree->SetBranchAddress("AliESDVertex",&vertex);
1150   Int_t entries = (Int_t)tree->GetEntries();
1151
1152   for(Int_t i=0; i<entries; i++) {
1153     total += 1;
1154     tree->GetEvent(i);
1155     nc = vertex->GetNContributors();
1156     if(nc<mincontr) continue;
1157     vertex->GetXYZ(vtx);
1158     vertex->GetSigmaXYZ(errvtx);
1159     //printf("%f %f %f %f %f %f\n",vtx[0],vtx[1],vtx[2],errvtx[0],errvtx[1],errvtx[2]);
1160
1161     avx += vtx[0];
1162     avy += vtx[1];
1163     avz += vtx[2];
1164     sum += 1.;
1165     wgtavx += vtx[0]/errvtx[0]/errvtx[0];
1166     wgtavy += vtx[1]/errvtx[1]/errvtx[1];
1167     wgtavz += vtx[2]/errvtx[2]/errvtx[2];
1168     sumwgtx += 1./errvtx[0]/errvtx[0];
1169     sumwgty += 1./errvtx[1]/errvtx[1];
1170     sumwgtz += 1./errvtx[2]/errvtx[2];
1171
1172     hx->Fill(1e4*vtx[0]);
1173     hy->Fill(1e4*vtx[1]);
1174     hz->Fill(vtx[2]);
1175
1176     vertex = 0;
1177   }
1178
1179   avx /= sum;
1180   avy /= sum;
1181   avz /= sum;
1182   wgtavx /= sumwgtx;
1183   wgtavy /= sumwgty;
1184   wgtavz /= sumwgtz;
1185   ewgtavx = 1./TMath::Sqrt(sumwgtx);
1186   ewgtavy = 1./TMath::Sqrt(sumwgty);
1187   ewgtavz = 1./TMath::Sqrt(sumwgtz);
1188
1189
1190   for(Int_t i=0; i<entries; i++) {
1191     tree->GetEvent(i);
1192     nc = vertex->GetNContributors();
1193     if(nc<mincontr) continue;
1194     vertex->GetXYZ(vtx);
1195     vertex->GetSigmaXYZ(errvtx);
1196
1197     varx += (vtx[0]-avx)*(vtx[0]-avx);
1198     vary += (vtx[1]-avy)*(vtx[1]-avy);
1199     varz += (vtx[2]-avz)*(vtx[2]-avz);
1200     covxy += (vtx[0]-avx)*(vtx[1]-avy);
1201     covxz += (vtx[0]-avx)*(vtx[2]-avz);
1202     covyz += (vtx[1]-avy)*(vtx[2]-avz);
1203
1204     vertex = 0;
1205   }
1206
1207   varx /= sum;
1208   vary /= sum;
1209   varz /= sum;
1210   covxy /= sum;
1211   covxz /= sum;
1212   covyz /= sum;
1213   rmsx = TMath::Sqrt(varx);
1214   rmsy = TMath::Sqrt(vary);
1215   rmsz = TMath::Sqrt(varz);
1216   eavx = rmsx/TMath::Sqrt(sum);
1217   eavy = rmsy/TMath::Sqrt(sum);
1218   eavz = rmsz/TMath::Sqrt(sum);
1219
1220
1221   printf("\n\nNumber of events: Total %d, Used %d\n",total,sum);
1222   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);
1223   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);
1224   printf("RMS:\n x = %f mum\n y = %f mum\n z = %f mum\n",1e4*rmsx,1e4*rmsy,1e4*rmsz);
1225
1226   TCanvas *c = new TCanvas("c","c",0,0,1000,500);
1227   c->Divide(3,1);
1228   c->cd(1);
1229   hx->Draw();
1230   TF1 *gx = new TF1("gx","gaus",-1000,1000);
1231   gx->SetLineColor(2);
1232   hx->Fit(gx);
1233   TF1 *gxx = (TF1*)gx->Clone("gxx");
1234   gxx->FixParameter(2,50.);
1235   gxx->SetLineStyle(2);
1236   gxx->Draw("same");
1237   c->cd(2);
1238   hy->Draw();
1239   TF1 *gy = new TF1("gy","gaus",-1000,1000);
1240   gy->SetLineColor(2);
1241   hy->Fit(gy);
1242   TF1 *gyy = (TF1*)gy->Clone("gyy");
1243   gyy->FixParameter(2,50.);
1244   gyy->SetLineStyle(2);
1245   gyy->Draw("same");
1246   c->cd(3);
1247   hz->Draw();
1248   TF1 *gz = new TF1("gz","gaus",-10,10);
1249   gz->SetLineColor(2);
1250   hz->Fit(gz);
1251   TF1 *gzz = (TF1*)gz->Clone("gzz");
1252   gzz->FixParameter(2,5.3);
1253   gzz->SetLineStyle(2);
1254   gzz->Draw("same");
1255
1256   vtx[0] = wgtavx;
1257   vtx[1] = wgtavy;
1258   vtx[2] = wgtavz;
1259   covvtx[0] = varx;
1260   covvtx[1] = covxy;
1261   covvtx[2] = vary;
1262   covvtx[3] = covxz;
1263   covvtx[4] = covyz;
1264   covvtx[5] = 5.3;
1265   AliESDVertex *vtxmean = new AliESDVertex(vtx,covvtx,0.,(Int_t)sum);
1266   vtxmean->SetTitle("Vertex from weighted average");
1267   TFile *out = new TFile("AliESDVertexMean.root","recreate");
1268   vtxmean->Write("vtxmean");
1269   out->Close();
1270
1271   return;
1272 }
1273 //---------------------------------------------------------------------------