ITS/AliITSFindPrimaryVertex.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TROOT.h>
   3 #include <TClassTable.h>
   4 #include <TStopwatch.h>
   5 #include <TFile.h>
   6 #include <TTree.h>
   7 #include <Riostream.h>
   8 #include <AliRun.h>
   9 #include <AliHeader.h>
  10 #include <AliGenEventHeader.h>
  11 #include <AliITSVertexerIons.h>
  12 #include <AliRunLoader.h>
  13 #include <AliITSVertexerIons.h>
  14 #include <AliITSLoader.h>
  15 #include <AliGenHijingEventHeader.h>
  16 #include <unistd.h>
  17
  18 #endif
  19
  20 void AliITSFindPrimaryVertex(Int_t evNumber1=0,Int_t NumbofEv=1, const char *filename="galice.root") {
  21
  22   Int_t evNumber2 = evNumber1+NumbofEv;
  23
  24   if (gClassTable->GetID("AliRun") < 0) {
  25     gROOT->Macro("loadlibs.C");
  26   } else {
  27     if(gAlice){
  28       delete AliRunLoader::Instance();
  29       delete gAlice;
  30       gAlice=0;
  31     }
  32   }
  33
  34   AliRunLoader* rl = AliRunLoader::Open("galice.root");
  35   if (rl == 0x0){
  36     ::Error("AliITSFindPrimaryVertex.C","Can not open session RL=NULL");
  37     return;
  38   }
  39
  40   // Open output file for vertices (default name: ITS.Vertex.root
  41   // and Create vertexer
  42
  43   AliITSVertexerIons *vertexer = new AliITSVertexerIons();
  44   vertexer->Init("default");
  45   //vertexer->SetDebug(1);
  46
  47   AliESDVertex *V;
  48   //   Loop over events
  49
  50   AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  51   itsloader->LoadRecPoints("read");
  52
  53   for (int nev=evNumber1; nev< evNumber2; nev++) {
  54     cout<<"=============================================================\n";
  55     cout<<" Processing event "<<nev<<endl;
  56     cout<<"=============================================================\n";
  57     rl->GetEvent(nev);
  58     AliHeader *header = rl->GetHeader();
  59     AliGenEventHeader* genEventHeader = header->GenEventHeader();
  60     TArrayF primaryVertex(3);
  61     genEventHeader->PrimaryVertex(primaryVertex);
  62
  63     AliGenHijingEventHeader* hijingHeader = (AliGenHijingEventHeader*)  genEventHeader;
  64     Float_t b = hijingHeader->ImpactParameter();
  65     cout << "Impact parameter = " << b << " fm" << endl;
  66
  67     TStopwatch timer;
  68     timer.Start();
  69
  70     TTree* cltree = itsloader->TreeR();
  71     V=vertexer->FindVertexForCurrentEvent(cltree);
  72
  73     TVector3 vtrue(primaryVertex[0],primaryVertex[1],primaryVertex[2]);
  74     TVector3 vfound(V->GetXv(),V->GetYv(),V->GetZv());
  75     TVector3 dif=vtrue-vfound;
  76     cout << "True vertex coordinates (cm) = " << vtrue.X() << " " << vtrue.Y() << " " << vtrue.Z() << endl;
  77     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;
  78
  79     timer.Stop();
  80     timer.Print();
  81
  82     vertexer->WriteCurrentVertex();
  83   }
  84
  85
  86 }
  87