Revised algorithm extended to very peripheral heavy-ion collisions (from G. Lo Re)
authorbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 11 Mar 2004 16:44:05 +0000 (16:44 +0000)
committerbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 11 Mar 2004 16:44:05 +0000 (16:44 +0000)
ITS/AliITSFindPrimaryVertex.C
ITS/AliITSVertexerIons.cxx
ITS/AliITSVertexerIons.h

index 095fcd1..388d0ea 100644 (file)
@@ -8,7 +8,9 @@
 #include <AliGenEventHeader.h>
 #include <AliITSVertexerIons.h>
 #include <AliRunLoader.h>
+#include <AliITSVertexerIons.h>
 #include <AliITSLoader.h>
+#include <unistd.h>
 
 #endif
 
@@ -25,48 +27,47 @@ void AliITSFindPrimaryVertex(Int_t evNumber1=0,Int_t NumbofEv=1, const char *fil
       gAlice=0;
     }
   }
-
   
   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;
-  }
-
+  
   // Open output file for vertices (default name: ITS.Vertex.root 
   // and Create vertexer
+
   AliITSVertexerIons *vertexer = new AliITSVertexerIons("default");
+  //vertexer->SetDebug(1);
+  
   AliESDVertex *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);
-  
+    
+    AliGenHijingEventHeader* hijingHeader = (AliGenHijingEventHeader*)  genEventHeader;
+    Float_t b = hijingHeader->ImpactParameter();   
+    cout << "Impact parameter = " << b << " fm" << endl;
+
     TStopwatch timer;
     timer.Start();
 
     V=vertexer->FindVertexForCurrentEvent(nev);
+
+    TVector3 vtrue(primaryVertex[0],primaryVertex[1],primaryVertex[2]);
+    TVector3 vfound(V->GetXv(),V->GetYv(),V->GetZv());
+    TVector3 dif=vtrue-vfound;
+    cout << "True vertex coordinates (cm) = " << vtrue.X() << " " << vtrue.Y() << " " << vtrue.Z() << endl;
+    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;
+    
     if(V){
       Double_t pos[3];
       for(Int_t kk=0;kk<3;kk++)pos[kk]=(Double_t)primaryVertex[kk];
@@ -74,17 +75,7 @@ void AliITSFindPrimaryVertex(Int_t evNumber1=0,Int_t NumbofEv=1, const char *fil
     }
     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;
-                
+    
     vertexer->WriteCurrentVertex(); 
   } 
   
index 0f95dcf..012c334 100644 (file)
 //                                                                  //
 //                                                                  //
 // Written by Giuseppe Lo Re and Francesco Riggi                    //
-//
 // Giuseppe.Lore@ct.infn.it                                         //
 // Franco.Riggi@ct.infn.it                                          //
 //                                                                  //
-// Release date: May 2001                                           //
+// Release date: Mar 2004                                           //
 //                                                                  //
 //                                                                  //       
 //////////////////////////////////////////////////////////////////////
@@ -44,12 +43,14 @@ ClassImp(AliITSVertexerIons)
 
 
 
-  //______________________________________________________________________
+//______________________________________________________________________
   AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
   // Default Constructor
 
   fITS = 0;
   SetNpThreshold();
+  SetMaxDeltaPhi();
+  SetMaxDeltaZ();
 }
 
 //______________________________________________________________________
@@ -58,7 +59,8 @@ AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
   
   fITS = 0;
   SetNpThreshold();
-}
+  SetMaxDeltaPhi();
+  SetMaxDeltaZ();}
 
 
 //______________________________________________________________________
@@ -69,14 +71,12 @@ AliITSVertexerIons::~AliITSVertexerIons() {
 
 //______________________________________________________________________
 AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
-  // Defines the AliESDVertex for the current event
+
   fCurrentVertex = 0;
-  Double_t position[3];
-  Double_t resolution[3];
-  Double_t snr[3];
+
   AliRunLoader *rl = AliRunLoader::GetRunLoader();
   AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
-  if(!fITS) { 
+  if(!fITS) {
     fITS =(AliITS *)gAlice->GetDetector("ITS");
     if(!fITS) {
       Error("FindVertexForCurrentEvent","AliITS object was not found");
@@ -84,105 +84,16 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
     }
   }
   fITS->SetTreeAddress();
-  AliITSgeom *g2 = fITS->GetITSgeom(); 
+  AliITSgeom *g2 = fITS->GetITSgeom();
   TClonesArray  *recpoints = fITS->RecPoints();
-  TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
-  TBranch *branch;
-  if(fUseV2Clusters){
-    branch = itsloader->TreeR()->GetBranch("Clusters");
-    branch->SetAddress(&clusters);
-  }
-  else {
-    branch = itsloader->TreeR()->GetBranch("ITSRecPoints");
-  }
-
   AliITSRecPoint *pnt;
+  TTree *tr =  itsloader->TreeR();
 
-  Int_t nopoints1 = 0;                                         
-  Int_t nopoints2 = 0;
-  Double_t vzero[3];
-  Double_t aspar[2];
-     
-  //------------ Rough Vertex evaluation ---------------------------------   
-
-  Int_t i,npoints,ipoint,j,k,max,binmax;
-  Double_t zCentroid;
+  Int_t npoints=0;
+  Int_t nopoints1=40000;
+  Int_t nopoints2=40000;
   Float_t l[3], p[3];
 
-  Double_t mxpiu = 0;
-  Double_t mxmeno = 0;
-  Double_t mypiu = 0;
-  Double_t mymeno = 0;
-   
-  TH1F *hITSz1         = new TH1F("hITSz1","",100,-14.35,14.35);
-  TTree *tr =  itsloader->TreeR(); 
-  for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) 
-    {
-      fITS->ResetRecPoints();
-      tr->GetEvent(i); 
-      if(fUseV2Clusters){
-       Clusters2RecPoints(clusters,i,recpoints);
-      }
-      npoints = recpoints->GetEntries();
-      for (ipoint=0;ipoint<npoints;ipoint++) {
-       pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
-       l[0]=pnt->GetX();
-       l[1]=0;
-       l[2]=pnt->GetZ();
-       g2->LtoG(i, l, p);
-       if(i<80 && TMath::Abs(p[2])<14.35) {         
-         if(p[0]>0) mxpiu++; 
-         if(p[0]<0) mxmeno++; 
-         if(p[1]>0) mypiu++; 
-         if(p[1]<0) mymeno++; 
-         hITSz1->Fill(p[2]);       
-       }
-       if(i>=80 && TMath::Abs(p[2])<14.35) nopoints2++;
-      }       
-    }  
-  
-  nopoints1 = (Int_t)(hITSz1->GetEntries()); 
-  if (nopoints1 == 0) {
-    delete hITSz1;
-    return fCurrentVertex;
-  }         
-  aspar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
-  aspar[1] = (mypiu-mymeno)/(mypiu+mymeno); 
-   
-  vzero[0] = 5.24441*aspar[0]; 
-  vzero[1] = 5.24441*aspar[1]; 
-
-  zCentroid=  TMath::Abs(hITSz1->GetMean()); 
-  vzero[2] = -5.31040e-02+1.42198*zCentroid+7.44718e-01*TMath::Power(zCentroid,2)
-    -5.73426e-01*TMath::Power(zCentroid,3)+2.01500e-01*TMath::Power(zCentroid,4)          
-    -3.34118e-02*TMath::Power(zCentroid,5)+2.20816e-03*TMath::Power(zCentroid,6);
-   
-  if(hITSz1->GetMean()<0) vzero[2] = -vzero[2];
-   
-  //cout << "\nXvzero: " << vzero[0] << " cm" << "";
-  //cout << "\nYvzero: " << vzero[1] << " cm" << "";
-  //cout << "\nZvzero: " << vzero[2] << " cm" << "\n";
-
-  delete hITSz1;
-
-  Double_t dphi,r,deltaZ,a,b;
-  Int_t np1=0;
-  Int_t np2=0;
-  Int_t niter=0;
-   
-  Double_t deltaPhiZ = 0.08;
-   
-  Float_t mag[3];
-  Float_t origin[3];
-  for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
-  gAlice->Field()->Field(origin,mag);
-
-  deltaPhiZ = deltaPhiZ*TMath::Abs(mag[2])/2;
-  Double_t deltaPhiXY = 1.0;   
-
-  //cout << "\ndeltaPhiZ: " << deltaPhiZ << " deg" << "\n";   
-  //cout << "deltaPhiXY: " << deltaPhiXY << " deg" << "\n";   
-   
   Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
   z1=new Double_t[nopoints1];
   z2=new Double_t[nopoints2];
@@ -194,276 +105,158 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
   phi2=new Double_t[nopoints2];
   r1=new Double_t[nopoints1];
   r2=new Double_t[nopoints2];
-        
-  deltaZ = 4.91617e-01+2.67567e-02*vzero[2]+1.49626e-02*TMath::Power(vzero[2],2); 
-  Float_t multFactorZ = 28000./(Float_t)nopoints1;
-  Int_t nbin=(Int_t)((deltaZ/0.005)/multFactorZ);
-  Int_t nbinxy=250;
-  Int_t *vectorBinZ,*vectorBinXY;
-  vectorBinZ=new Int_t[nbin];
-  vectorBinXY=new Int_t[nbinxy];
-  Float_t f1= 0;
-  Float_t f2= 0;
-  Double_t sigma,averagebg;
-     
-  TH1D *hITSZv         = new TH1D("hITSZv","",nbin,vzero[2]-deltaZ,vzero[2]+deltaZ);
-  TH1D *hITSXv         = new TH1D("hITSXv","",nbinxy,-3,3);
-  TH1D *hITSYv         = new TH1D("hITSYv","",nbinxy,-3,3);
-
-  //cout << "deltaZeta: " << deltaZ << " cm" << "\n";   
-   
-   
- start:   
-
-  hITSZv->Add(hITSZv,-1);
-  hITSXv->Add(hITSXv,-1);
-  hITSYv->Add(hITSYv,-1);
 
-  np1=np2=0;
-
-
-  for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) 
-    {
-      fITS->ResetRecPoints(); 
-      itsloader->TreeR()->GetEvent(i);
-     if(fUseV2Clusters){
-       Clusters2RecPoints(clusters,i,recpoints);
+  Double_t mxpiu = 0;
+  Double_t mxmeno = 0;
+  Double_t mypiu = 0;
+  Double_t mymeno = 0;
+  Double_t r=0;
+
+  Int_t np1=0, np2=0;
+  for(Int_t i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) {
+    fITS->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[1]=0;
+      l[2]=pnt->GetZ();
+      g2->LtoG(i, l, p);
+      r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
+
+      if(i<80 && TMath::Abs(p[2])<14.35)  {
+               y1[np1]=p[1];
+               x1[np1]=p[0];
+               z1[np1]=p[2];
+               if(p[0]>0) mxpiu++;
+               if(p[0]<0) mxmeno++;
+               if(p[1]>0) mypiu++;
+               if(p[1]<0) mymeno++;
+               np1++;
       }
-      npoints = recpoints->GetEntries();
-      for (ipoint=0;ipoint<npoints;ipoint++) {
-               
-       pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
-       l[0]=pnt->GetX();
-       l[1]=0;
-       l[2]=pnt->GetZ();
-       g2->LtoG(i, l, p);
-              
-       if(i<80 && TMath::Abs(p[2])<14.35) {
-         p[0]=p[0]-vzero[0];
-         p[1]=p[1]-vzero[1];
-         r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
-         y1[np1]=p[1];
-         x1[np1]=p[0];
-         z1[np1]=p[2]; 
-         r1[np1]=r;
-         phi1[np1]=PhiFunc(p);
-         np1++;
-       }
-
-       if(i>=80 &&  TMath::Abs(p[2])<14.35) { 
-         p[0]=p[0]-vzero[0];
-         p[1]=p[1]-vzero[1];
-         r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
-         y2[np2]=p[1];
-         x2[np2]=p[0];
-         z2[np2]=p[2];
-         r2[np2]=r;
-         phi2[np2]=PhiFunc(p);
-         np2++;
-       }
-                        
+      if(i>=80 && TMath::Abs(p[2])<14.35) {
+               y2[np2]=p[1];
+               x2[np2]=p[0];
+               z2[np2]=p[2];
+               np2++;
       }
     }
+  }
 
-  //------------------ Correlation between rec points ----------------------
-    
-  
   if(np1<fNpThreshold) {
-    /*    cout << "Warning: AliITSVertexerIons finder is not reliable for low multiplicity events. May be you have to use AliITSVertexerPPZ.\n";
-         position[0]=vzero[0];
-         position[1]=vzero[1];
-         position[2]=vzero[2];
-         resolution[0]=-123;
-         resolution[1]=-123;
-         resolution[2]=-123;
-         snr[0]=-123;
-         snr[1]=-123;
-         snr[2]=-123;
-         Char_t name[30];
-         sprintf(name,"Vertex");
-         fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
-         return fCurrentVertex;*/
-
     Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerPPZ with default parameters...\n");
     Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
-    AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("null");
+    AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("default");
     fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
     delete dovert;
     return fCurrentVertex;
   }
 
-
-  for(j=0; j<(np2)-1; j++) {
-    for(k=0; k<(np1)-1; k++) { 
-      dphi=TMath::Abs(phi1[k]-phi2[j]);
-      if(dphi>180) dphi = 360-dphi;
-      if(dphi<deltaPhiZ && TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-vzero[2])
-        <deltaZ) hITSZv->Fill(z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1));        
-    }
-  }
-      
-  //cout << "\nNumber of used pairs: \n";
-
-  a = vzero[2]-deltaZ; 
-  b = vzero[2]+deltaZ; 
-  max=(Int_t) hITSZv->GetMaximum();
-  binmax=hITSZv->GetMaximumBin();
-  sigma=0;
-  for(i=0;i<nbin;i++) vectorBinZ[i]=(Int_t)hITSZv->GetBinContent(i);
-  for(i=0;i<10;i++) f1=f1+vectorBinZ[i]/10;
-  for(i=nbin-10;i<nbin;i++) f2=f2+vectorBinZ[i]/10;
-  averagebg=(f1+f2)/2;
-  for(i=0;i<nbin;i++) {
-    if(vectorBinZ[i]-averagebg>(max-averagebg)*0.4 && 
-       vectorBinZ[i]-averagebg<(max-averagebg)*0.7) {
-      sigma=hITSZv->GetBinCenter(binmax)-hITSZv->GetBinCenter(i);
-      sigma=TMath::Abs(sigma);
-      if(sigma==0) sigma=0.05;
-    }
-  }
-  
-  
-  TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);  
-  
-  fz->SetParameter(0,max);
-  if(niter==0) {Double_t temp = hITSZv->GetBinCenter(binmax); vzero[2]=temp;}
-  fz->SetParameter(1,vzero[2]);
-  fz->SetParameter(2,sigma); 
-  fz->SetParameter(3,averagebg);  
-  fz->SetParLimits(2,0,999);
-  fz->SetParLimits(3,0,999);
-  
-  hITSZv->Fit("fz","RQ0");
-  
-  snr[2] = fz->GetParameter(0)/fz->GetParameter(3);
-  if(snr[2]<0.) { 
-    Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for z!!!\n");
-    Error("FindVertexForCurrentEvent","The algorithm cannot find the z vertex position.\n");
-    //    exit(123456789);
+  if(!np1 || !np2) {
+    Error("FindVertexForCurrentEvent","No points in the pixer layers");
     return fCurrentVertex;
+  }  
+  if(fDebug>0)   cout << "N. points layer 1 and 2 : " << np1 << " " << np2 << endl;
+
+  Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
+  Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
+
+  Double_t x0 = 2.86*asparx;
+  Double_t y0 = 2.86*aspary;
+
+  if(fDebug>0)   cout << "Rough xy vertex = " << x0 << " " << y0 << endl;  
+
+  for(Int_t i=0;i<np1;i++) {
+    x1[i]-=x0;
+    y1[i]-=y0;
+    PhiFunc(x1[i],y1[i],phi1[i]);
+    r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]);
   }
-  else {
-    position[2] = fz->GetParameter(1);
-    if(position[2]<0) 
-      {
-       position[2]=position[2]-TMath::Abs(position[2])*1.11/10000;
-      }
-    else
-      {
-       position[2]=position[2]+TMath::Abs(position[2])*1.11/10000;
-      }
+  for(Int_t i=0;i<np2;i++) {
+    x2[i]-=x0;
+    y2[i]-=y0;
+    PhiFunc(x2[i],y2[i],phi2[i]);
+    r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
   }
-  resolution[2] = fz->GetParError(1);
-  
-  for(j=0; j<(np2)-1; j++) {
-    for(k=0; k<(np1)-1; k++) { 
-      dphi=TMath::Abs(phi1[k]-phi2[j]);
+   
+  Int_t nbinxy=400;
+  Int_t nbinz=400;
+  TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5);
+  TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5);
+  TH1F *hzv=new TH1F("hzv","",nbinz,-15,15);
+
+  Double_t dphi;
+  for(Int_t j=0; j<np1; j++) {
+    for(Int_t k=0; k<np2; k++) {
+      dphi=TMath::Abs(phi2[k]-phi1[j]);
       if(dphi>180) dphi = 360-dphi;
-      if(dphi>deltaPhiXY) continue;
-      if(TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-position[2])
-        <4*resolution[2]) {
-       hITSXv->Fill(vzero[0]+(x2[j]-((x2[j]-x1[k])/(y2[j]-y1[k]))*y2[j]));
-       hITSYv->Fill(vzero[1]+(y2[j]-((y2[j]-y1[k])/(x2[j]-x1[k]))*x2[j]));
-      }
-    }
-  }
-  
-  TF1 *fx = new TF1("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",vzero[0]-0.5,vzero[0]+0.5);  
-  
-  max=(Int_t) hITSXv->GetMaximum();
-  binmax=hITSXv->GetMaximumBin();
-  sigma=0;
-  f1=f2=0;
-  for(i=0;i<nbinxy;i++) vectorBinXY[i]=(Int_t)hITSXv->GetBinContent(i);
-  for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10;
-  for(i=nbinxy-10;i<nbinxy;i++) f2=f2+vectorBinXY[i]/10;
-  averagebg=(f1+f2)/2;
-  for(i=0;i<nbinxy;i++) {
-    if(vectorBinXY[i]-averagebg>(max-averagebg)*0.4 && 
-       vectorBinXY[i]-averagebg<(max-averagebg)*0.7) {
-      sigma=hITSXv->GetBinCenter(binmax)-hITSXv->GetBinCenter(i);
-      sigma=TMath::Abs(sigma);
-      if(sigma==0) sigma=0.05;
-    }
+      if(dphi>fMaxDeltaPhi) continue; 
+      hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
+     }
   }
+
+  // Fitting ...
+  Double_t max;
+  Int_t bin_max;
+  Double_t max_center;
   
+  max = hzv->GetMaximum();
+  bin_max=hzv->GetMaximumBin();
+  max_center=hzv->GetBinCenter(bin_max);
+  Double_t dxy=1.5;
   
-  fx->SetParameter(0,max);
-  fx->SetParameter(1,vzero[0]);
-  fx->SetParameter(2,sigma); 
-  fx->SetParameter(3,averagebg);
-  
-  hITSXv->Fit("fx","RQ0"); 
-  
-  snr[0] = fx->GetParameter(0)/fx->GetParameter(3);
-  
-  if(snr[0]<0.) { 
-    Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for x!\n");
-    Error("FindVertexForCurrentEvent","The algorithm cannot find the x vertex position\n");
-    //   exit(123456789);  
+  TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",max_center-dxy,max_center+dxy);  
+  fz->SetParameter(0,max);
+  fz->SetParameter(1,max_center);
+  fz->SetParameter(2,0.1);
+  fz->SetLineColor(kRed);
+  hzv->Fit("fz","RQ0");
+
+  if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) {
+    Warning("FindVertexForCurrentEvent","AliITSVertexerIonsGeneral finder is not reliable for events with x and y vertex coordinates too far from the origin of the reference system. Their values will be set to 0.\n");
+    Double_t position[3]={0,0,fz->GetParameter(1)};
+    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);
+    sprintf(name,"Vertex");
+    fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
     return fCurrentVertex;
   }
-  else {
-    position[0]=fx->GetParameter(1);
-    resolution[0]=fx->GetParError(1);
-  }
-  
-  TF1 *fy = new TF1("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",vzero[1]-0.5,vzero[1]+0.5);  
   
-  max=(Int_t) hITSYv->GetMaximum();
-  binmax=hITSYv->GetMaximumBin();
-  sigma=0;
-  f1=f2=0;
-  for(i=0;i<nbinxy;i++) vectorBinXY[i]=(Int_t)hITSYv->GetBinContent(i);
-  for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10;
-  for(i=nbinxy-10;i<nbinxy;i++) f2=f2+vectorBinXY[i]/10;
-  averagebg=(f1+f2)/2;
-  for(i=0;i<nbinxy;i++) {
-    if(vectorBinXY[i]-averagebg>(max-averagebg)*0.4 && 
-       vectorBinXY[i]-averagebg<(max-averagebg)*0.7) {
-      sigma=hITSYv->GetBinCenter(binmax)-hITSYv->GetBinCenter(i);
-      sigma=TMath::Abs(sigma);
-      if(sigma==0) sigma=0.05;
+  for(Int_t j=0; j<np1; j++) {
+    for(Int_t k=0; k<np2; k++) {
+      if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue;
+      if(y2[k]==y1[j]) continue;
+      hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k]));
+      if(x2[k]==x1[j]) continue;
+      hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k]));
     }
   }
+    
+  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);
+  fx->SetParameter(2,0.1);
+  fx->SetLineColor(kRed);
+  hxv->Fit("fx","RQW0"); 
   
-  fy->SetParameter(0,max);
-  fy->SetParameter(1,vzero[1]);
-  fy->SetParameter(2,sigma); 
-  fy->SetParameter(3,averagebg);
-  
-  hITSYv->Fit("fy","RQ0"); 
-  
-  snr[1] = fy->GetParameter(0)/fy->GetParameter(3);
-  if(snr[1]<0.) {   
-    Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for y!\n");
-    Error("FindVertexForCurrentEvent","The algorithm cannot find the y vertex position.\n");
-    //    exit(123456789); 
-    return fCurrentVertex;
-  }
-  else {
-    position[1]=fy->GetParameter(1);
-    resolution[1]=fy->GetParError(1);
-  }
+  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);
+  fy->SetParameter(2,0.1);
+  fy->SetLineColor(kRed);
+  hyv->Fit("fy","RQW0");
   
-    
-  /*cout << "iter = " << niter << " -> x = " << position[0] << endl;
-    cout << "iter = " << niter << " -> y = " << position[1] << endl;
-    cout << "iter = " << niter << " -> z = " << position[2] << endl;
-    cout << "---\n";*/
-
-  niter++;
-   
-  vzero[0] = position[0];
-  vzero[1] = position[1];
-  vzero[2] = position[2];
-
-  if(niter<3) goto start;  // number of iterations
-   
-
-  cout<<"Number of iterations: "<<niter<<endl;
-  cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl; 
   delete [] z1;
   delete [] z2;
   delete [] y1;
@@ -474,14 +267,15 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
   delete [] r2;
   delete [] phi1;
   delete [] phi2;
-  delete [] vectorBinZ;
-  delete [] vectorBinXY;
-   
-  delete hITSZv;
-  delete hITSXv;
-  delete hITSYv;
+  delete hxv;
+  delete hyv;
+  delete hzv;
+  
+  Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)};
+  Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)};
+  Double_t snr[3]={0,0,0};
+  
   Char_t name[30];
-  //  sprintf(name,"Vertex_%d",evnumber);
   if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
   sprintf(name,"Vertex");
   fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
@@ -489,14 +283,11 @@ AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
 }
 
 //______________________________________________________________________
-Double_t AliITSVertexerIons::PhiFunc(Float_t p[]) {
-  Double_t phi=0;
-
-  if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578);
-  if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
-  if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
-  if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360;
-  return phi;
+void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
+  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;
+  if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;;
 }
 
 //______________________________________________________________________
@@ -528,3 +319,19 @@ void AliITSVertexerIons::PrintStatus() const {
   cout<<"\n Last event to be processed "<<fLastEvent<<endl;
   if(fCurrentVertex)fCurrentVertex->PrintStatus();
 }
+
+//________________________________________________________
+Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
+  Int_t max_content=0;
+  Int_t max_bin=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;}
+  }
+  Double_t max=h->GetBinCenter(max_bin);
+  return max;
+}
+
+
index 5f4c0d5..67d7fcc 100644 (file)
@@ -2,9 +2,10 @@
 #define ALIITSVERTEXERIONS_H
 
 #include <AliITSVertexer.h>
+#include <TH1F.h>
 
 //////////////////////////////////////////////////////////////////////
-// AliITSVertexerIons is a class for full 3D primary vertex         //
+// AliITSVertexerIons  is a class for full 3D primary vertex        //
 // finding optimized for Ion-Ion interactions                       //
 //                                                                  // 
 //                                                                  //
 //                                                                  //
 // Written by Giuseppe Lo Re and Francesco Riggi                    //
 // Giuseppe.Lore@ct.infn.it                                         //
-//                                                                  //
 // Franco.Riggi@ct.infn.it                                          //
 //                                                                  //
-// Release date: May 2001                                           //
+// Release date: Mar 2004                                           //
 //                                                                  //
 //                                                                  //       
 //////////////////////////////////////////////////////////////////////
@@ -30,16 +30,23 @@ class AliITSVertexerIons : public AliITSVertexer {
   virtual ~AliITSVertexerIons(); // destructor
   virtual AliESDVertex* FindVertexForCurrentEvent(Int_t event);
   virtual void FindVertices();
-  virtual Double_t PhiFunc(Float_t p[]);
+  virtual void PhiFunc(Double_t &x,Double_t &y,Double_t &phi);
   virtual void PrintStatus() const;
   Int_t GetNpThreshold() const {return fNpThreshold;}
-  void SetNpThreshold(Int_t t = 4000){fNpThreshold = t;}
+  void SetNpThreshold(Int_t t = 500){fNpThreshold = t;}
+  Double_t GetMaxDeltaPhi() const {return fMaxDeltaPhi;}
+  void SetMaxDeltaPhi(Double_t dphi=0.45) {fMaxDeltaPhi=dphi;}
+  Double_t GetMaxDeltaZ() const {return fMaxDeltaPhi;}
+  void SetMaxDeltaZ(Double_t dz=0.15) {fMaxDeltaZ=dz;}
+  Double_t FindMaxAround(Double_t point, TH1F *h, Double_t distance);
 
  protected:
   AliITS *fITS;            //! pointer to the AliITS object
   Int_t fNpThreshold;      // minimum number of rec points for vertexing
+  Double_t fMaxDeltaPhi;
+  Double_t fMaxDeltaZ;
 
-  ClassDef(AliITSVertexerIons,2);
+  ClassDef(AliITSVertexerIons,3);
 };
 
 #endif