From 504de69bcd1df138d148e077f995e02ef7f8aaf2 Mon Sep 17 00:00:00 2001 From: barbera Date: Tue, 27 Mar 2001 14:51:48 +0000 Subject: [PATCH] New class for vertex finding --- ITS/AliITSVertex.cxx | 311 +++++++++++++++++++++++++++++++++++++++++++ ITS/AliITSVertex.h | 47 +++++++ 2 files changed, 358 insertions(+) create mode 100644 ITS/AliITSVertex.cxx create mode 100644 ITS/AliITSVertex.h diff --git a/ITS/AliITSVertex.cxx b/ITS/AliITSVertex.cxx new file mode 100644 index 00000000000..846c685daaf --- /dev/null +++ b/ITS/AliITSVertex.cxx @@ -0,0 +1,311 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "AliRun.h" +#include "AliITS.h" +#include "AliITSgeom.h" +#include "AliITSRecPoint.h" +#include "AliGenerator.h" +#include "AliMagF.h" + +#include +#include +#include +#include +#include +#include + + +ClassImp(AliITSVertex) + +//////////////////////////////////////////////////////////////////////// +// AliITSVertex is a class for primary vertex finding. +// +// Version: 0 +// Written by Giuseppe Lo Re and Francesco Riggi +// Giuseppe.Lore@ct.infn.it +// Franco.Riggi@ct.infn.it +// Marh 2001 +// +// +/////////////////////////////////////////////////////////////////////// + + +//______________________________________________________________________________ + +AliITSVertex::AliITSVertex() { + + fPosition = new Double_t[3]; + fResolution = new Double_t[3]; + fSNR = new Double_t[3]; + + AliITS* aliits =(AliITS *)gAlice->GetDetector("ITS"); + AliITSgeom *g2 = ((AliITS*)aliits)->GetITSgeom(); + + TClonesArray *recpoints = aliits->RecPoints(); + AliITSRecPoint *pnt; + TRandom rnd; + + Int_t NoPoints1 = 0; + Int_t NoPoints2 = 0; + Double_t Vzero[3]; + Double_t AsPar[2]; + +//------------Rough Vertex------------------------------------------------------ + + Int_t i,npoints,ipoint,j,k; + Double_t ZCentroid; + Float_t l[3], p[3]; + + TH1F *hITSx1pos = new TH1F("hITSx1pos","",100,0.,4.2); + TH1F *hITSx1neg = new TH1F("hITSx1neg","",100,-4.2,0.); + TH1F *hITSy1pos = new TH1F("hITSy1pos","",100,0.,4.2); + TH1F *hITSy1neg = new TH1F("hITSy1neg","",100,-4.2,0.); + TH1F *hITSz1 = new TH1F("hITSz1","",100,-14.35,14.35); + + for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) + { + aliits->ResetRecPoints(); + gAlice->TreeR()->GetEvent(i); + + npoints = recpoints->GetEntries(); + for (ipoint=0;ipointUncheckedAt(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) hITSx1pos->Fill(p[0]); + if(p[0]<0) hITSx1neg->Fill(p[0]); + if(p[1]>0) hITSy1pos->Fill(p[1]); + if(p[1]<0) hITSy1neg->Fill(p[1]); + hITSz1->Fill(p[2]); + } + if(i>=80 && TMath::Abs(p[2])<14.35) (NoPoints2)++; + } + } + + NoPoints1 = (Int_t)(hITSz1->GetEntries()); + + Double_t mxpiu = hITSx1pos->GetEntries(); + Double_t mxmeno = hITSx1neg->GetEntries(); + Double_t mypiu = hITSy1pos->GetEntries(); + Double_t mymeno = hITSy1neg->GetEntries(); + + AsPar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno); + AsPar[1] = (mypiu-mymeno)/(mypiu+mymeno); + + Vzero[0] = 5.24434*AsPar[0]; + Vzero[1] = 5.24434*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 << "\nCentroid and RMS of hIITSz1: \n"; + //cout << hITSz1->GetMean() << " " << hITSz1->GetRMS() <GetStartSPD();i<=g2->GetLastSPD();i++) + { + aliits->ResetRecPoints(); + gAlice->TreeR()->GetEvent(i); + npoints = recpoints->GetEntries(); + for (ipoint=0;ipoint2) continue; + pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint); + l[0]=pnt->GetX(); + l[1]=0; + l[2]=pnt->GetZ(); + g2->LtoG(i, l, p); + 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)); + + if(i<80 && TMath::Abs(p[2])<14.35) { + Z1[np1]=p[2]; + Y1[np1]=p[1]; + X1[np1]=p[0]; + r1[np1]=r; + phi1[np1]=PhiFunc(p); + np1++; + } + + if(i>=80 && TMath::Abs(p[2])<14.35) { + Z2[np2]=p[2]; + Y2[np2]=p[1]; + X2[np2]=p[0]; + r2[np2]=r; + phi2[np2]=PhiFunc(p); + np2++; + } + + } + } + +//-------------------Correlations----------------------------------------------- + + //cout << "\nPoints number on the first and second layer: "<Field()->Field(origin,B); + + DeltaPhi = DeltaPhi*B[2]/2; + + cout << "\nDeltaPhi: " << DeltaPhi << " deg" << "\n"; + + DeltaZ = + 0.2+2.91617e-01+2.67567e-02*Vzero[2]+1.49626e-02*TMath::Power(Vzero[2],2); + + cout << "DeltaZeta: " << DeltaZ << " cm" << "\n"; + + TH1F *hITSZv = new TH1F("hITSZv","",DeltaZ/0.005,Vzero[2]-DeltaZ,Vzero[2]+DeltaZ); + + for(j=0; j<(np2)-1; j++) { + for(k=0; k<(np1)-1; k++) { + if(TMath::Abs(Z1[k]-Z2[j])>13.) continue; + dphi=TMath::Abs(phi1[k]-phi2[j]); + if(dphi>180) dphi = 360-dphi; + if(dphiFill(Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1)); + } + } + } + + + + cout << "\nNumber of used pairs: \n"; + cout << hITSZv->GetEntries() << '\n' << '\n'; + cout << "\nCentroid and RMS of hIITSZv: \n"; + cout << hITSZv->GetMean() << " " << hITSZv->GetRMS() << "\n"<< "\n"; + + delete [] Z1; + delete [] Z2; + delete [] Y1; + delete [] Y2; + delete [] X1; + delete [] X2; + delete [] r1; + delete [] r2; + delete [] phi1; + delete [] phi2; + + Double_t a = Vzero[2]-DeltaZ; + Double_t b = Vzero[2]+DeltaZ; + + TF1 *f4 = new TF1 ("f4","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b); + f4->SetParameter(0,700.); + f4->SetParameter(1,Vzero[2]); + f4->SetParameter(2,0.05); // !! + f4->SetParameter(3,50.); + + hITSZv->Fit("f4","RME0"); + + delete hITSZv; + + fSNR[2] = f4->GetParameter(0)/f4->GetParameter(3); + if(fSNR[2]<1.5) { + cout << "\nSignal to noise ratio too small!!!" << endl; + cout << "The algorithm can't find the z vertex position." << endl; + exit(123456789); + } + else + { + fPosition[2] = (f4->GetParameter(1)); + if(fPosition[2]<0) + { + fPosition[2]=fPosition[2]-TMath::Abs(fPosition[2])*1.11/10000; + } + else + { + fPosition[2]=fPosition[2]+TMath::Abs(fPosition[2])*1.11/10000; + } + } + fResolution[2] = f4->GetParError(1); + + AliGenerator *gener = gAlice->Generator(); + Float_t Vx,Vy,Vz; + gener->GetOrigin(Vx,Vy,Vz); + + fPosition[0]=(Double_t)Vx; + fPosition[1]=(Double_t)Vy; + + fResolution[0] = 0; + fResolution[1] = 0; + + fSNR[0] = 0; + fSNR[1] = 0; + +} + + +//______________________________________________________________________________ + +AliITSVertex::~AliITSVertex() { + delete [] fPosition; + delete [] fResolution; + delete [] fSNR; +} + +//______________________________________________________________________________ + + +Double_t AliITSVertex::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; +} +// +//______________________________________________________________________________ + diff --git a/ITS/AliITSVertex.h b/ITS/AliITSVertex.h new file mode 100644 index 00000000000..c7ea24ae918 --- /dev/null +++ b/ITS/AliITSVertex.h @@ -0,0 +1,47 @@ +#ifndef ALIITSVERTEX_H +#define ALIITSVERTEX_H + +class TTree; +class TFile; +class AliITSgeom; +class AliITSRecPoint; +class TH1F; +class TF1; +class TClonesArray; +class TObject; +class AliGenerator; + +class AliITSVertex : public TObject { + + public: + + AliITSVertex(); + ~AliITSVertex(); + Double_t PhiFunc(Float_t p[]); + +// At present this class determines vertex position, resolution and signal +// to noise ratio only for the z coordinate. For x and y coordinates it +// gives the values x and y setted in Config.C with resolution and signal +// to noise ratio values = 0. +// The cases of beam off-set and magnetic field = 0.4 T are included. + + Double_t GetZv() {return (Double_t)fPosition[2];} + Double_t GetZRes() {return fResolution[2];} + Double_t GetZSNR() {return fSNR[2];} + Double_t GetYv() {return (Double_t)fPosition[1];} + Double_t GetYRes() {return fResolution[1];} + Double_t GetYSNR() {return fSNR[1];} + Double_t GetXv() {return (Double_t)fPosition[0];} + Double_t GetXRes() {return fResolution[0];} + Double_t GetXSNR() {return fSNR[0];} + + private: + + Double_t *fPosition; + Double_t *fResolution; + Double_t *fSNR; + +ClassDef(AliITSVertex,1) // Class for Vertex finder +}; + +#endif -- 2.43.5