Changes related to the extraction of the V0 finder into a separate class (A. Dainese...

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

diff --git a/ITS/AliITSFindPrimaryVertex.C b/ITS/AliITSFindPrimaryVertex.C
index c3c72e16a8f7832594088d05f8085417156c671c..5e343a2fd182b31d40842a716f8058cac46bac1b 100644 (file)
--- a/ITS/AliITSFindPrimaryVertex.C
+++ b/ITS/AliITSFindPrimaryVertex.C
@@ -1,5 +1,7 @@
 #if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TROOT.h>
 #include <TClassTable.h>
+#include <TStopwatch.h>
 #include <TFile.h>
 #include <TTree.h>
 #include <Riostream.h>
@@ -8,7 +10,10 @@
 #include <AliGenEventHeader.h>
 #include <AliITSVertexerIons.h>
 #include <AliRunLoader.h>
+#include <AliITSVertexerIons.h>
 #include <AliITSLoader.h>
+#include <AliGenHijingEventHeader.h>
+#include <unistd.h>
 
 #endif
 
@@ -20,71 +25,60 @@ void AliITSFindPrimaryVertex(Int_t evNumber1=0,Int_t NumbofEv=1, const char *fil
     gROOT->Macro("loadlibs.C");
   } else {
     if(gAlice){
-      delete gAlice->GetRunLoader();
+      delete AliRunLoader::Instance();
       delete gAlice;
       gAlice=0;
     }
   }
-
   
   AliRunLoader* rl = AliRunLoader::Open("galice.root");
   if (rl == 0x0){
     ::Error("AliITSFindPrimaryVertex.C","Can not open session RL=NULL");
     return;
   }
- 
-  Int_t retval = rl->LoadHeader();
-  if (retval){
-    cerr<<"AliITSFindPrimaryVertex.C : LoadHeader returned error"<<endl;
-    return;
-  }
-  retval = rl->LoadKinematics();
-  if (retval){
-    cerr<<"AliITSFindPrimaryVertex.C : LoadKinematics returned error"<<endl;
-    return;
-  }
-
+  
   // Open output file for vertices (default name: ITS.Vertex.root 
   // and Create vertexer
-  AliITSVertexerIons *vertexer = new AliITSVertexerIons("default");
-  AliITSVertex *V;
+
+  AliITSVertexerIons *vertexer = new AliITSVertexerIons();
+  vertexer->Init("default");
+  //vertexer->SetDebug(1);
+  
+  AliESDVertex *V;
   //   Loop over events 
-  //
- 
+   
+  AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
+  itsloader->LoadRecPoints("read");
+
   for (int nev=evNumber1; nev< evNumber2; nev++) {
     cout<<"=============================================================\n";
     cout<<" Processing event "<<nev<<endl;
     cout<<"=============================================================\n";
     rl->GetEvent(nev);
-    cout << "nev         " << nev <<endl;
-    // The true Z coord. is fetched for comparison
     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;
+
     TStopwatch timer;
     timer.Start();
 
-    V=vertexer->FindVertexForCurrentEvent(nev);
-    if(V){
-      Double_t pos[3];
-      for(Int_t kk=0;kk<3;kk++)pos[kk]=(Double_t)primaryVertex[kk];
-      V->SetTruePos(pos);
-    }
+    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();
-
-    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;
-                
+    
     vertexer->WriteCurrentVertex(); 
   }