Some changes for new coordinate system

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

diff --git a/ITS/AliITSFindPrimaryVertex.C b/ITS/AliITSFindPrimaryVertex.C
index fd8b8229e879a700f8d4bed37c2e936bc8a069b2..39576aa0d8582ab14bc386f29459915bda113de2 100644 (file)
--- a/ITS/AliITSFindPrimaryVertex.C
+++ b/ITS/AliITSFindPrimaryVertex.C
@@ -1,66 +1,93 @@
-#include "iostream.h"
-#include "TMath.h"
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TClassTable.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 <AliITSLoader.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 gAlice->GetRunLoader();
+      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;
-  
   
-     TStopwatch timer;
-     timer.Start();
-
-     AliITSVertex *V = new AliITSVertex();
-         V->Exec();
+  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;
+  }
 
-     timer.Stop();
-     timer.Print();
+  // Open output file for vertices (default name: ITS.Vertex.root 
+  // and Create vertexer
+  AliITSVertexerIons *vertexer = new AliITSVertexerIons("default");
+  AliITSVertex *V;
+  //   Loop over events 
+  //
+ 
+  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);
+  
+    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;
+    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);
+    }
+    timer.Stop();
+    timer.Print();
+    if(!V)continue;
+    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;
+    vertexer->WriteCurrentVertex(); 
   } 
   
-  file->Close(); 
+
 }