ITS/AliITSFindPrimaryVertex.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TClassTable.h>
   3 #include <TFile.h>
   4 #include <TTree.h>
   5 #include <Riostream.h>
   6 #include <AliRun.h>
   7 #include <AliHeader.h>
   8 #include <AliGenEventHeader.h>
   9 #include <AliITSVertexerIons.h>
  10 #include <AliRunLoader.h>
  11 #include <AliITSLoader.h>
  12
  13 #endif
  14
  15 void AliITSFindPrimaryVertex(Int_t evNumber1=0,Int_t NumbofEv=1, const char *filename="galice.root") {
  16
  17   Int_t evNumber2 = evNumber1+NumbofEv;
  18
  19   if (gClassTable->GetID("AliRun") < 0) {
  20     gROOT->Macro("loadlibs.C");
  21   } else {
  22     if(gAlice){
  23       delete gAlice->GetRunLoader();
  24       delete gAlice;
  25       gAlice=0;
  26     }
  27   }
  28
  29
  30   AliRunLoader* rl = AliRunLoader::Open("galice.root");
  31   if (rl == 0x0){
  32     ::Error("AliITSFindPrimaryVertex.C","Can not open session RL=NULL");
  33     return;
  34   }
  35
  36   Int_t retval = rl->LoadHeader();
  37   if (retval){
  38     cerr<<"AliITSFindPrimaryVertex.C : LoadHeader returned error"<<endl;
  39     return;
  40   }
  41   retval = rl->LoadKinematics();
  42   if (retval){
  43     cerr<<"AliITSFindPrimaryVertex.C : LoadKinematics returned error"<<endl;
  44     return;
  45   }
  46
  47   // Open output file for vertices (default name: ITS.Vertex.root
  48   // and Create vertexer
  49   AliITSVertexerIons *vertexer = new AliITSVertexerIons("default");
  50   AliITSVertex *V;
  51   //   Loop over events
  52   //
  53
  54   for (int nev=evNumber1; nev< evNumber2; nev++) {
  55     cout<<"=============================================================\n";
  56     cout<<" Processing event "<<nev<<endl;
  57     cout<<"=============================================================\n";
  58     rl->GetEvent(nev);
  59     cout << "nev         " << nev <<endl;
  60     // The true Z coord. is fetched for comparison
  61     AliHeader *header = rl->GetHeader();
  62     AliGenEventHeader* genEventHeader = header->GenEventHeader();
  63     TArrayF primaryVertex(3);
  64     genEventHeader->PrimaryVertex(primaryVertex);
  65
  66     TStopwatch timer;
  67     timer.Start();
  68
  69     V=vertexer->FindVertexForCurrentEvent(nev);
  70     if(V){
  71       Double_t pos[3];
  72       for(Int_t kk=0;kk<3;kk++)pos[kk]=(Double_t)primaryVertex[kk];
  73       V->SetTruePos(pos);
  74     }
  75     timer.Stop();
  76     timer.Print();
  77     if(!V)continue;
  78     cout << endl << "Xv = " << V->GetXv() << " cm" << endl;
  79     cout << "X resolution = " << V->GetXRes()*10000 << " microns"  << endl;
  80     cout << "Signal/Noise for X = " << V->GetXSNR() << endl;
  81     cout << endl << "Yv = " << V->GetYv() << " cm"  << endl;
  82     cout << "Y resolution = " << V->GetYRes()*10000 << " microns"  << endl;
  83     cout << "Signal/Noise for Y = " << V->GetYSNR() << endl;
  84     cout << endl << "Zv = " << V->GetZv() << " cm" << endl;
  85     cout << "Z Resolution = " << V->GetZRes()*10000 << " microns" << endl;
  86     cout << "Signal/Noise for Z = " << V->GetZSNR() <<endl;
  87
  88     vertexer->WriteCurrentVertex();
  89   }
  90
  91
  92 }
  93