modified histo ranges

[u/mrichter/AliRoot.git] / ITS / AliITSFindPrimaryVertex.C

diff --git a/ITS/AliITSFindPrimaryVertex.C b/ITS/AliITSFindPrimaryVertex.C
index fd8b8229e879a700f8d4bed37c2e936bc8a069b2..5e343a2fd182b31d40842a716f8058cac46bac1b 100644 (file)
--- a/ITS/AliITSFindPrimaryVertex.C
+++ b/ITS/AliITSFindPrimaryVertex.C
@@ -1,66 +1,87 @@
-#include "iostream.h"
-#include "TMath.h"
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TROOT.h>
+#include <TClassTable.h>
+#include <TStopwatch.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <Riostream.h>
+#include <AliRun.h>
+#include <AliHeader.h>
+#include <AliGenEventHeader.h>
+#include <AliITSVertexerIons.h>
+#include <AliRunLoader.h>
+#include <AliITSVertexerIons.h>
+#include <AliITSLoader.h>
+#include <AliGenHijingEventHeader.h>
+#include <unistd.h>
 
-void AliITSFindPrimaryVertex(Int_t evNumber1=0,Int_t evNumber2=0, const char *filename="galice.root") {
+#endif
 
-  
-  ///////////////// Dynamically link some shared libs ////////////////////////////////
+void AliITSFindPrimaryVertex(Int_t evNumber1=0,Int_t NumbofEv=1, const char *filename="galice.root") {
+
+  Int_t evNumber2 = evNumber1+NumbofEv;
   
   if (gClassTable->GetID("AliRun") < 0) {
-    gROOT->LoadMacro("loadlibs.C");
-    loadlibs();
+    gROOT->Macro("loadlibs.C");
   } else {
-    delete gAlice;
-    gAlice=0;
+    if(gAlice){
+      delete AliRunLoader::Instance();
+      delete gAlice;
+      gAlice=0;
+    }
   }
-
-// Connect the Root Galice file containing Geometry, Kine and Hits
-   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
-   if (!file) file = new TFile(filename,"UPDATE");
-
-
-// Get AliRun object from file or create it if not on file
-   if (!gAlice) {
-      gAlice = (AliRun*)file->Get("gAlice");
-      if (gAlice) printf("AliRun object found on file\n");
-      if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
-   }
-
-//   Loop over events 
-//
-   Int_t Nh=0;
-   Int_t Nh1=0;
-   for (int nev=0; nev<= evNumber2; nev++) {
-     gAlice->SetEvent(nev);
-        Int_t nparticles = gAlice->GetEvent(nev);
-     cout << "nev         " << nev <<endl;
-     cout << "nparticles  " << nparticles <<endl;
-     if (nev < evNumber1) continue;
-     if (nparticles <= 0) return;
   
+  AliRunLoader* rl = AliRunLoader::Open("galice.root");
+  if (rl == 0x0){
+    ::Error("AliITSFindPrimaryVertex.C","Can not open session RL=NULL");
+    return;
+  }
+  
+  // Open output file for vertices (default name: ITS.Vertex.root 
+  // and Create vertexer
+
+  AliITSVertexerIons *vertexer = new AliITSVertexerIons();
+  vertexer->Init("default");
+  //vertexer->SetDebug(1);
   
-     TStopwatch timer;
-     timer.Start();
+  AliESDVertex *V;
+  //   Loop over events 
+   
+  AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
+  itsloader->LoadRecPoints("read");
 
-     AliITSVertex *V = new AliITSVertex();
-         V->Exec();
+  for (int nev=evNumber1; nev< evNumber2; nev++) {
+    cout<<"=============================================================\n";
+    cout<<" Processing event "<<nev<<endl;
+    cout<<"=============================================================\n";
+    rl->GetEvent(nev);
+    AliHeader *header = rl->GetHeader();
+    AliGenEventHeader* genEventHeader = header->GenEventHeader();
+    TArrayF primaryVertex(3);
+    genEventHeader->PrimaryVertex(primaryVertex);
+    
+    AliGenHijingEventHeader* hijingHeader = (AliGenHijingEventHeader*)  genEventHeader;
+    Float_t b = hijingHeader->ImpactParameter();   
+    cout << "Impact parameter = " << b << " fm" << endl;
 
-     timer.Stop();
-     timer.Print();
+    TStopwatch timer;
+    timer.Start();
 
-     cout << endl << "Xv = " << V->GetXv() << " cm" << endl;
-     cout << "X resolution = " << V->GetXRes()*10000 << " microns"  << endl;
-     cout << "Signal/Noise for X = " << V->GetXSNR() << endl;
-     cout << endl << "Yv = " << V->GetYv() << " cm"  << endl;
-     cout << "Y resolution = " << V->GetYRes()*10000 << " microns"  << endl;
-     cout << "Signal/Noise for Y = " << V->GetYSNR() << endl;
-     cout << endl << "Zv = " << V->GetZv() << " cm" << endl;
-     cout << "Z Resolution = " << V->GetZRes()*10000 << " microns" << endl;
-     cout << "Signal/Noise for Z = " << V->GetZSNR() <<endl;
-                
-     delete V;
+    TTree* cltree = itsloader->TreeR();
+    V=vertexer->FindVertexForCurrentEvent(cltree);
+
+    TVector3 vtrue(primaryVertex[0],primaryVertex[1],primaryVertex[2]);
+    TVector3 vfound(V->GetXv(),V->GetYv(),V->GetZv());
+    TVector3 dif=vtrue-vfound;
+    cout << "True vertex coordinates (cm) = " << vtrue.X() << " " << vtrue.Y() << " " << vtrue.Z() << endl;
+    cout << "Found vertex coordinates  (cm) = " << vfound.X() << " " << vfound.Y() << " " << vfound.Z() << endl;    cout << "Difference true - found (cm) = " << dif.Mag() << " " << dif.X() << " " << dif.Y() << " " << dif.Z() << endl;
+    
+    timer.Stop();
+    timer.Print();
+    
+    vertexer->WriteCurrentVertex(); 
   } 
   
-  file->Close(); 
+
 }