]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ITS/AliITSVertexerIons.cxx
Fix of parsing bug related to the reading of the calib header. Added consistency...
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerIons.cxx
index 012c334e07724d2a76873117133c859013b7dfee..eef9f19ac0a93bf7ceba306f2680256529835a11 100644 (file)
@@ -1,28 +1,16 @@
-#include <AliITSVertexerIons.h>
-#include "stdlib.h"
-#include <TMath.h>
-#include <TRandom.h>
-#include <TObjArray.h>
-#include <TROOT.h>
-#include <TFile.h>
-#include <TTree.h>
-#include <Riostream.h>
-
-#include "AliRun.h"
-#include "AliITS.h"
-#include "AliITSgeom.h"
+#include "AliRunLoader.h"
+#include "AliITSDetTypeRec.h"
+#include "AliITSVertexerIons.h"
+#include "AliITSVertexerZ.h"
+#include "AliESDVertex.h"
+#include "AliITSgeomTGeo.h"
 #include "AliITSLoader.h"
 #include "AliITSRecPoint.h"
-#include "AliGenerator.h"
-#include "AliMagF.h"
-
+#include "AliLog.h"
+#include <TTree.h>
 #include <TH1.h>
 #include <TF1.h>
-#include <TCanvas.h>
-#include <AliESDVertex.h>
-#include <TObjArray.h>
-#include <TObject.h>
-#include <AliITSVertexerPPZ.h>
+
 //////////////////////////////////////////////////////////////////////
 // AliITSVertexerIons is a class for full 3D primary vertex         //
 // finding optimized for Ion-Ion interactions                       //
 //                                                                  //       
 //////////////////////////////////////////////////////////////////////
 
-ClassImp(AliITSVertexerIons)
-
-
 
+ClassImp(AliITSVertexerIons)
+  
 //______________________________________________________________________
-  AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
+  AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer(),
+fNpThreshold(0),
+fMaxDeltaPhi(0),
+fMaxDeltaZ(0){
   // Default Constructor
 
-  fITS = 0;
+  //fITS = 0;
   SetNpThreshold();
   SetMaxDeltaPhi();
   SetMaxDeltaZ();
 }
 
 //______________________________________________________________________
-AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
+AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn),
+fNpThreshold(0),
+fMaxDeltaPhi(0),
+fMaxDeltaZ(0) {
   // Standard constructor
   
-  fITS = 0;
+  //fITS = 0;
   SetNpThreshold();
   SetMaxDeltaPhi();
-  SetMaxDeltaZ();}
+  SetMaxDeltaZ();
+}
+/*
+//______________________________________________________________________
+AliITSVertexerIons::AliITSVertexerIons(const AliITSVertexerIons &source):AliITSVertexer(source) {
+  // Copy constructor
+  // Copies are not allowed. The method is protected to avoid misuse.
+  Error("AliITSVertexerIons","Copy constructor not allowed\n");
+}
+*/
+//_________________________________________________________________________
+//AliITSVertexerIons& AliITSVertexerIons::operator=(const AliITSVertexerIons &/*source*/) {
+  // Assignment operator
+  // Assignment is not allowed. The method is protected to avoid misuse.
+  //Error("= operator","Assignment operator not allowed\n");
+  //return *this;
+//}
 
 
 //______________________________________________________________________
 AliITSVertexerIons::~AliITSVertexerIons() {
   // Default Destructor
-  fITS = 0;
+  //fITS = 0;
 }
 
 //______________________________________________________________________
-AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
+AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){ 
+// Defines the AliESDVertex for the current event
 
   fCurrentVertex = 0;
 
   AliRunLoader *rl = AliRunLoader::GetRunLoader();
-  AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
-  if(!fITS) {
-    fITS =(AliITS *)gAlice->GetDetector("ITS");
-    if(!fITS) {
-      Error("FindVertexForCurrentEvent","AliITS object was not found");
-      return fCurrentVertex;
-    }
-  }
-  fITS->SetTreeAddress();
-  AliITSgeom *g2 = fITS->GetITSgeom();
-  TClonesArray  *recpoints = fITS->RecPoints();
-  AliITSRecPoint *pnt;
+  AliITSLoader* itsloader = (AliITSLoader*)rl->GetLoader("ITSLoader");
+  /*
+  TDirectory * olddir = gDirectory;
+  rl->CdGAFile();
+  AliITSgeom* g2  = (AliITSgeom*)gDirectory->Get("AliITSgeom");
+  olddir->cd(); 
+  */
+
   TTree *tr =  itsloader->TreeR();
+  AliITSDetTypeRec detTypeRec;
+
+  detTypeRec.SetTreeAddressR(tr);
+
+  TClonesArray  *recpoints = detTypeRec.RecPoints();
+  AliITSRecPoint *pnt;
+
 
   Int_t npoints=0;
   Int_t nopoints1=40000;
   Int_t nopoints2=40000;
-  Float_t l[3], p[3];
+  Float_t p[3];
 
   Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
   z1=new Double_t[nopoints1];
@@ -113,16 +126,19 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
   Double_t r=0;
 
   Int_t np1=0, np2=0;
-  for(Int_t i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) {
-    fITS->ResetRecPoints();
+  for(Int_t i=AliITSgeomTGeo::GetModuleIndex(1,1,1);i<=AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;i++) {
+    detTypeRec.ResetRecPoints();
     tr->GetEvent(i);
     npoints = recpoints->GetEntries();
     for (Int_t ipoint=0;ipoint<npoints;ipoint++) {
       pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
-      l[0]=pnt->GetX();
+      /*
+      l[0]=pnt->GetDetLocalX();
       l[1]=0;
-      l[2]=pnt->GetZ();
+      l[2]=pnt->GetDetLocalZ();
       g2->LtoG(i, l, p);
+      */
+      pnt->GetGlobalXYZ(p);
       r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
 
       if(i<80 && TMath::Abs(p[2])<14.35)  {
@@ -145,9 +161,9 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
   }
 
   if(np1<fNpThreshold) {
-    Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerPPZ with default parameters...\n");
+    Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerZ with default parameters...\n");
     Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
-    AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("default");
+    AliITSVertexerZ *dovert = new AliITSVertexerZ("default");
     fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
     delete dovert;
     return fCurrentVertex;
@@ -157,7 +173,7 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
     Error("FindVertexForCurrentEvent","No points in the pixer layers");
     return fCurrentVertex;
   }  
-  if(fDebug>0)   cout << "N. points layer 1 and 2 : " << np1 << " " << np2 << endl;
+  AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2));
 
   Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
   Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
@@ -165,7 +181,7 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
   Double_t x0 = 2.86*asparx;
   Double_t y0 = 2.86*aspary;
 
-  if(fDebug>0)   cout << "Rough xy vertex = " << x0 << " " << y0 << endl;  
+  AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0));  
 
   for(Int_t i=0;i<np1;i++) {
     x1[i]-=x0;
@@ -198,17 +214,17 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
 
   // Fitting ...
   Double_t max;
-  Int_t bin_max;
-  Double_t max_center;
+  Int_t binMax;
+  Double_t maxcenter;
   
   max = hzv->GetMaximum();
-  bin_max=hzv->GetMaximumBin();
-  max_center=hzv->GetBinCenter(bin_max);
+  binMax=hzv->GetMaximumBin();
+  maxcenter=hzv->GetBinCenter(binMax);
   Double_t dxy=1.5;
   
-  TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",max_center-dxy,max_center+dxy);  
+  TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy);  
   fz->SetParameter(0,max);
-  fz->SetParameter(1,max_center);
+  fz->SetParameter(1,maxcenter);
   fz->SetParameter(2,0.1);
   fz->SetLineColor(kRed);
   hzv->Fit("fz","RQ0");
@@ -219,7 +235,7 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
     Double_t resolution[3]={0,0,0};
     Double_t snr[3]={0,0,0};
     Char_t name[30];
-    if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
+    AliDebug(1,Form("Vertex found for event %d",evnumber));
     sprintf(name,"Vertex");
     fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
     return fCurrentVertex;
@@ -238,21 +254,21 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
   TF1 *fx = new TF1 ("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",x0-dxy,x0+dxy);  
   fx->SetParameter(0,100);
   Double_t dist=0.3;
-  Double_t x_approx=FindMaxAround(x0,hxv,dist);
-  if(fDebug>0) cout << "x_approx = " << x_approx << endl;
-  fx->SetParameter(1,x_approx);
-  Double_t dif_centroid=0.07;
-  fx->SetParLimits(1,x_approx-dif_centroid,x_approx+dif_centroid);
+  Double_t xapprox=FindMaxAround(x0,hxv,dist);
+  AliDebug(1,Form("xapprox = %f",xapprox));
+  fx->SetParameter(1,xapprox);
+  Double_t difcentroid=0.07;
+  fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid);
   fx->SetParameter(2,0.1);
   fx->SetLineColor(kRed);
   hxv->Fit("fx","RQW0"); 
   
   TF1 *fy = new TF1 ("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",y0-dxy,y0+dxy);  
   fy->SetParameter(0,100);
-  Double_t y_approx=FindMaxAround(y0,hyv,dist);
-  if(fDebug>0) cout << "y_approx = " << y_approx << endl;
-  fy->SetParameter(1,y_approx);
-  fy->SetParLimits(1,y_approx-dif_centroid,y_approx+dif_centroid);
+  Double_t yapprox=FindMaxAround(y0,hyv,dist);
+  AliDebug(1,Form("yapprox = %f",yapprox));
+  fy->SetParameter(1,yapprox);
+  fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid);
   fy->SetParameter(2,0.1);
   fy->SetLineColor(kRed);
   hyv->Fit("fy","RQW0");
@@ -276,14 +292,16 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
   Double_t snr[3]={0,0,0};
   
   Char_t name[30];
-  if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
+  AliDebug(1,Form("Vertex found for event %d",evnumber));
   sprintf(name,"Vertex");
   fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
+
   return fCurrentVertex;
 }
 
 //______________________________________________________________________
 void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
+  // Method for the computation of the phi angle given x and y. The result is in degrees. 
   if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578);
   if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
   if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
@@ -303,9 +321,7 @@ void AliITSVertexerIons::FindVertices(){
       WriteCurrentVertex();
     }
     else {
-      if(fDebug>0){
-       cout<<"Vertex not found for event "<<i<<endl;
-      }
+      AliDebug(1,Form("Vertex not found for event %d",i));
     }
   }
 }
@@ -314,7 +330,6 @@ void AliITSVertexerIons::FindVertices(){
 void AliITSVertexerIons::PrintStatus() const {
   // Print current status
   cout <<"=======================================================\n";
-  cout <<" Debug flag: "<<fDebug<<endl;
   cout<<"First event to be processed "<<fFirstEvent;
   cout<<"\n Last event to be processed "<<fLastEvent<<endl;
   if(fCurrentVertex)fCurrentVertex->PrintStatus();
@@ -322,15 +337,16 @@ void AliITSVertexerIons::PrintStatus() const {
 
 //________________________________________________________
 Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
-  Int_t max_content=0;
-  Int_t max_bin=0;
+  // It finds the maximum of h within point-distance and distance+point. 
+  Int_t maxcontent=0;
+  Int_t maxbin=0;
   for(Int_t i=0;i<h->GetNbinsX();i++) {
     Int_t content=(Int_t)h->GetBinContent(i);
     Double_t center=(Double_t)h->GetBinCenter(i);
-    if(fabs(center-point)>distance) continue;
-    if(content>max_content) {max_content=content;max_bin=i;}
+    if(TMath::Abs(center-point)>distance) continue;
+    if(content>maxcontent) {maxcontent=content;maxbin=i;}
   }
-  Double_t max=h->GetBinCenter(max_bin);
+  Double_t max=h->GetBinCenter(maxbin);
   return max;
 }