From c95d6417c3345e0344913fabc0741ace5cc9276f Mon Sep 17 00:00:00 2001
From: masera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Fri, 21 Feb 2003 16:02:22 +0000
Subject: [PATCH] Algorithm efficiency improved

---
 ITS/AliITSVertexerPPZ.cxx | 118 +++++++++++++++++++++++++++++++++++---
 ITS/AliITSVertexerPPZ.h   |  11 +++-
 2 files changed, 119 insertions(+), 10 deletions(-)

diff --git a/ITS/AliITSVertexerPPZ.cxx b/ITS/AliITSVertexerPPZ.cxx
index 6c676b816ec..32f0765a9b9 100644
--- a/ITS/AliITSVertexerPPZ.cxx
+++ b/ITS/AliITSVertexerPPZ.cxx
@@ -1,3 +1,17 @@
+/**************************************************************************
+ * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
 #include <AliITSVertexerPPZ.h>
 #include <Riostream.h>
 #include <TArrayF.h>
@@ -36,6 +50,7 @@ AliITSVertexerPPZ::AliITSVertexerPPZ():AliITSVertexer() {
   fZFound = 0;
   fZsig = 0.;
   fITS = 0;
+  SetWindow(0);
 }
 
 AliITSVertexerPPZ::AliITSVertexerPPZ(TFile *infile, TFile *outfile, Float_t x0, Float_t y0):AliITSVertexer(infile,outfile) {
@@ -49,6 +64,7 @@ AliITSVertexerPPZ::AliITSVertexerPPZ(TFile *infile, TFile *outfile, Float_t x0,
   fZFound = 0;
   fZsig = 0.;
   fITS = 0;
+  SetWindow();
 }
 
 //______________________________________________________________________
@@ -59,6 +75,7 @@ AliITSVertexerPPZ::~AliITSVertexerPPZ() {
 
 //________________________________________________________
 void  AliITSVertexerPPZ::EvalZ(TH1F *hist,Int_t sepa, Int_t ncoinc, TArrayF *zval) {
+  Float_t DeltaVal = hist->GetBinWidth(1)*fWindow; // max window in Z for searching
   fZFound=0;
   fZsig=0;
   Int_t N=0;
@@ -68,28 +85,55 @@ void  AliITSVertexerPPZ::EvalZ(TH1F *hist,Int_t sepa, Int_t ncoinc, TArrayF *zva
   Float_t curz = 0.;
   for(Int_t i=1;i<=sepa;i++){
     Float_t cont=hist->GetBinContent(i);
-    curz = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(i);
+    if(cont!=0)curz = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(i);
     totst+=cont;
     totst2+=cont*cont;
     N++;
     if(cont!=0)NbinNotZero++;
   }
   if(NbinNotZero==0){fZFound=-100; fZsig=-100; return;}
-  if(NbinNotZero==1){fZFound = curz; fZsig=-100; return;}
-  totst2=TMath::Sqrt(totst2/(N-1)-totst*totst/N/(N-1));
+  if(NbinNotZero==1){
+    fZFound = curz;
+    fZsig=0;
+    fCurrentVertex = new AliITSVertex(fZFound,fZsig,NbinNotZero);
+    return;
+  }
+  Float_t errsq = totst2/(N-1)-totst*totst/N/(N-1);
+  if(errsq>=0){
+  totst2=TMath::Sqrt(errsq);
+  }
+  else {
+    Error("EvalZ","Negative variance: %d - %12.7f - %12.7f",N,totst2,totst);
+    fZFound=-100; 
+    fZsig=-100;
+    return;
+  }
   totst /= N;
+  Float_t cut = totst+totst2*2.;
+  if(fDebug>1)cout<<"totst, totst2, cut: "<<totst<<", "<<totst2<<", "<<cut<<endl;
   Float_t val1=hist->GetBinLowEdge(sepa); 
   Float_t val2=hist->GetBinLowEdge(1);
+  Float_t valm = 0.;
+  Float_t zmax = 0.;
   for(Int_t i=1;i<=sepa;i++){
     Float_t cont=hist->GetBinContent(i);
-    if(cont>(totst+totst2*2.)){
+    if(cont>valm){
+      valm = cont;
+      zmax = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(1);
+    }
+    if(cont>cut){
       curz=hist->GetBinLowEdge(i);
       if(curz<val1)val1=curz;
       if(curz>val2)val2=curz;
     }
   }
   val2+=hist->GetBinWidth(1);
-  if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<endl;
+  if((val2-val1)>DeltaVal){
+    val1 = zmax-DeltaVal/2.;
+    val2 = zmax+DeltaVal/2.;
+    if(fDebug>0)cout<<"val1 and val2 recomputed\n";
+  }
+  if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<" "<<val2-val1<<endl;
   fZFound=0;
   fZsig=0;
   N=0;
@@ -97,13 +141,36 @@ void  AliITSVertexerPPZ::EvalZ(TH1F *hist,Int_t sepa, Int_t ncoinc, TArrayF *zva
     Float_t z=zval->At(i);
     if(z<val1)continue;
     if(z>val2)continue;
+   
     fZFound+=z;
     fZsig+=z*z;
+    /*   weights defined by the curvature
+    Float_t wei = 1./curv->At(i);
+    fZFound+=z*wei;
+    fZsig+=wei;
+    */
     N++;
   }
-  if(N<=1){fZFound=-110; fZsig=-110; return;}
-  fZsig=TMath::Sqrt((fZsig/(N-1)-fZFound*fZFound/N/(N-1))/N);
+  if(N<1){fZFound=-110; fZsig=-110; return;}
+  if(N==1){
+    fZsig = 0;
+    fCurrentVertex = new AliITSVertex(fZFound,fZsig,N);
+    return;
+  }
+  errsq = (fZsig/(N-1)-fZFound*fZFound/N/(N-1))/N;
+  if(errsq>=0.){
+    fZsig=TMath::Sqrt(errsq);
+  }
+  else {
+    Error("evalZ","Negative variance: %d - %12.7f %12.7f",N,fZsig,fZFound);
+    fZsig=0;
+  }
   fZFound=fZFound/N;
+  /* weights defined by the curvature
+  fZsig=1./fZsig;
+  fZFound*=fZsig;
+  fZsig = TMath::Sqrt(fZsig);
+  */
   fCurrentVertex = new AliITSVertex(fZFound,fZsig,N);
 }
 
@@ -174,6 +241,7 @@ AliITSVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
   if(firipixe>1){
     rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
     zave=zave/firipixe;
+    if(fDebug>1)cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++\n Number of firing pixels: "<<firipixe<<endl;
   }
   else {
     fZFound = -200;
@@ -186,9 +254,13 @@ AliITSVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
   Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
   zlim2=zlim1 + sepa/10.;
   TH1F *zvdis = new TH1F("z_ev","zv distr",sepa,zlim1,zlim2);
-  if(fDebug>0)cout<<"Z limits: "<<zlim1<<" "<<zlim2<<"; Bins= "<<sepa<<endl;
+  if(fDebug>0){
+    cout<<"Z limits: "<<zlim1<<" "<<zlim2<<"; Bins= "<<sepa<<endl;
+    cout<<"Bin width: "<<zvdis->GetBinWidth(1)*10000.<<" micron\n";
+  }
   Int_t sizarr=100;
   TArrayF *zval = new TArrayF(sizarr);
+  //  TArrayF *curv = new TArrayF(sizarr);
   Int_t ncoinc=0;
   for(Int_t module= fFirstL1; module<=fLastL1;module++){
     if(fDebug>0)cout<<"processing module   "<<module<<"                  \r";
@@ -228,10 +300,16 @@ AliITSVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
 	    if(diff<fDiffPhiMax ){
 	      zvdis->Fill(zr0);
 	      zval->AddAt(zr0,ncoinc);
+	      /* uncomment these lines to use curvature as a weight
+	      Float_t cu = Curv(0.,0.,gc[0],gc[1],gc2[0],gc[1]);
+	      curv->AddAt(cu,ncoinc);
+	      */
 	      ncoinc++;
 	      if(ncoinc==(sizarr-1)){
 		sizarr+=100;
 		zval->Set(sizarr);
+		//uncomment next line to use curvature as weight
+		//		curv->Set(sizarr);
 	      }
 	    }
 	  }
@@ -244,9 +322,11 @@ AliITSVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
   if(fDebug>0){
     cout<<endl<<"Number of coincidences = "<<ncoinc<<endl;
   }
+  //  EvalZ(zvdis,sepa,ncoinc,zval,curv);
   EvalZ(zvdis,sepa,ncoinc,zval);
   delete zvdis;
   delete zval;
+  //  delete curv;
   if(fCurrentVertex){
     char name[30];
     sprintf(name,"Vertex_%d",evnumber);
@@ -295,6 +375,7 @@ void AliITSVertexerPPZ::PrintStatus() const {
   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
   cout <<fLastL2<<endl;
   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
+  cout <<" Window for Z search: "<<fWindow<<endl;
   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
   cout <<" Debug flag: "<<fDebug<<endl;
   if(fInFile)cout <<" The input file is open\n";
@@ -304,3 +385,24 @@ void AliITSVertexerPPZ::PrintStatus() const {
   cout<<"First event to be processed "<<fFirstEvent;
   cout<<"\n Last event to be processed "<<fLastEvent<<endl;
 }
+
+//________________________________________________________
+Float_t AliITSVertexerPPZ::Curv(Double_t x1,Double_t y1,
+                   Double_t x2,Double_t y2,
+                   Double_t x3,Double_t y3) 
+{
+  //-----------------------------------------------------------------
+  // Initial approximation of the track curvature (Y. Belikov) squared
+  //-----------------------------------------------------------------
+  Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
+  Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
+                  (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
+  Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
+                  (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
+
+  Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
+
+  Float_t val = (Float_t)(( -xr*yr/TMath::Sqrt(xr*xr+yr*yr))
+			  *( -xr*yr/TMath::Sqrt(xr*xr+yr*yr)));
+  return val; 
+}
diff --git a/ITS/AliITSVertexerPPZ.h b/ITS/AliITSVertexerPPZ.h
index f6e31769e2a..c2ff716e3ac 100644
--- a/ITS/AliITSVertexerPPZ.h
+++ b/ITS/AliITSVertexerPPZ.h
@@ -1,5 +1,8 @@
 #ifndef ALIITSVERTEXERPPZ_H
 #define ALIITSVERTEXERPPZ_H
+/* Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
 
 #include <AliITSVertexer.h>
 
@@ -25,9 +28,12 @@ class AliITSVertexerPPZ : public AliITSVertexer {
   virtual Float_t GetZFound() const {return fZFound;}
   virtual Float_t GetZsig() const {return fZsig;}
   virtual void PrintStatus() const;
-  virtual void SetDiffPhiMax(Float_t pm = 0.01){fDiffPhiMax = pm;}
+  virtual void SetDiffPhiMax(Float_t pm = 0.05){fDiffPhiMax = pm;}
   virtual void SetFirstLayerModules(Int_t m1 = 0, Int_t m2 = 79){fFirstL1 = m1; fLastL1 = m2;}
   virtual void SetSecondLayerModules(Int_t m1 = 80, Int_t m2 = 239){fFirstL2 = m1; fLastL2 = m2;}
+  virtual void SetWindow(Float_t w=3.){fWindow = w;}
+  static Float_t Curv(Double_t x1,Double_t y1, Double_t x2,Double_t y2,
+		       Double_t x3,Double_t y3); 
 
  private:
   void EvalZ(TH1F *hist,Int_t sepa, Int_t ncoinc, TArrayF *zval);
@@ -43,9 +49,10 @@ class AliITSVertexerPPZ : public AliITSVertexer {
   AliITS *fITS;            //! pointer to the AliITS object
   Float_t fZFound;         //! found value for the current event
   Float_t fZsig;           //! RMS of Z
+  Float_t fWindow;         // window width for Z search in mm (3 mm by def.)
 
 
-  ClassDef(AliITSVertexerPPZ,1);
+  ClassDef(AliITSVertexerPPZ,2);
 };
 
 #endif
-- 
2.39.3