1 /**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 #include <AliITSVertexerZ.h>
16 #include <Riostream.h>
19 #include "AliITSLoader.h"
22 #include<TClonesArray.h>
26 #include <AliITSgeom.h>
27 #include <AliITSRecPoint.h>
29 /////////////////////////////////////////////////////////////////
30 // this class implements a fast method to determine
31 // the Z coordinate of the primary vertex
32 // for p-p collisions it seems to give comparable or better results
33 // with respect to what obtained with AliITSVertexerPPZ
34 // It can be used successfully with Pb-Pb collisions
35 ////////////////////////////////////////////////////////////////
37 ClassImp(AliITSVertexerZ)
41 //______________________________________________________________________
42 AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer() {
43 // Default constructor
47 SetFirstLayerModules(0);
48 SetSecondLayerModules(0);
56 SetBinWidthCoarse(0.);
61 //______________________________________________________________________
62 AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
63 // Standard Constructor
67 SetFirstLayerModules();
68 SetSecondLayerModules();
82 //______________________________________________________________________
83 AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr) {
85 // Copies are not allowed. The method is protected to avoid misuse.
86 Error("AliITSVertexerZ","Copy constructor not allowed\n");
89 //______________________________________________________________________
90 AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ& /* vtxr */){
91 // Assignment operator
92 // Assignment is not allowed. The method is protected to avoid misuse.
93 Error("= operator","Assignment operator not allowed\n");
98 //______________________________________________________________________
99 AliITSVertexerZ::~AliITSVertexerZ() {
100 // Default Destructor
102 if(fZCombc)delete fZCombc;
103 if(fZCombf)delete fZCombf;
106 //______________________________________________________________________
107 AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
108 // Defines the AliESDVertex for the current event
111 AliRunLoader *rl =AliRunLoader::GetRunLoader();
112 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
113 itsLoader->ReloadRecPoints();
115 rl->GetEvent(evnumber);
118 fITS = (AliITS*)gAlice->GetModule("ITS");
120 Error("FindVertexForCurrentEvent","AliITS object was not found");
121 return fCurrentVertex;
125 fITS->SetTreeAddress();
127 AliITSgeom *geom = fITS->GetITSgeom();
129 TTree *tR = itsLoader->TreeR();
130 TClonesArray *itsRec = 0;
131 Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
132 Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
133 Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
134 Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
136 itsRec = fITS->RecPoints();
137 TBranch *branch = tR->GetBranch("ITSRecPoints");
139 Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
140 Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
141 if(fZCombc)delete fZCombc;
142 fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
143 if(fZCombf)delete fZCombf;
144 fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
148 for(Int_t module= fFirstL1; module<=fLastL1;module++){
149 if(module%4==0 || module%4==3)continue;
150 branch->GetEvent(module);
151 nrpL1+= itsRec->GetEntries();
152 fITS->ResetRecPoints();
154 for(Int_t module= fFirstL2; module<=fLastL2;module++){
155 branch->GetEvent(module);
156 nrpL2+= itsRec->GetEntries();
157 fITS->ResetRecPoints();
159 cout<<"nrpL1 = "<<nrpL1<<"; nrpL2= "<<nrpL2<<endl;
160 if(nrpL1 == 0 || nrpL2 == 0){
161 cout<<"achi, achi\n";
163 return fCurrentVertex;
165 Float_t *xc1 = new Float_t [nrpL1];
166 Float_t *yc1 = new Float_t [nrpL1];
167 Float_t *zc1 = new Float_t [nrpL1];
168 Float_t *phi1 = new Float_t [nrpL1];
169 Float_t *xc2 = new Float_t [nrpL2];
170 Float_t *yc2 = new Float_t [nrpL2];
171 Float_t *zc2 = new Float_t [nrpL2];
172 Float_t *phi2 = new Float_t [nrpL2];
174 for(Int_t module= fFirstL1; module<=fLastL1;module++){
175 if(module%4==0 || module%4==3)continue;
176 branch->GetEvent(module);
177 Int_t nrecp1 = itsRec->GetEntries();
178 for(Int_t j=0;j<nrecp1;j++){
179 AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
182 geom->LtoG(module,lc,gc);
188 phi1[ind] = TMath::ATan2(gc[1],gc[0]);
189 if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
192 fITS->ResetRecPoints();
195 for(Int_t module= fFirstL2; module<=fLastL2;module++){
196 branch->GetEvent(module);
197 Int_t nrecp2 = itsRec->GetEntries();
198 for(Int_t j=0;j<nrecp2;j++){
199 AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
202 geom->LtoG(module,lc,gc);
208 phi2[ind] = TMath::ATan2(gc[1],gc[0]);
209 if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
212 fITS->ResetRecPoints();
214 for(Int_t i=0;i<nrpL1;i++){
215 Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]);
216 for(Int_t j=0;j<nrpL2;j++){
217 Float_t diff = TMath::Abs(phi2[j]-phi1[i]);
218 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
219 if(diff<fDiffPhiMax){
220 Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]);
221 Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1);
235 if(fZCombc->GetEntries()==0){
236 Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
238 return fCurrentVertex;
241 cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
243 Int_t bi = fZCombc->GetMaximumBin();
244 Float_t centre = fZCombc->GetBinCenter(bi);
245 Int_t n1 = static_cast<Int_t>((centre-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
246 Int_t n2 = static_cast<Int_t>((centre+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
254 for(Int_t n=n1;n<=n2;n++){
255 fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
256 num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
257 fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
263 Float_t radi = fZsig/(num-1)-fZFound*fZFound/num/(num-1);
264 if(radi>0.)fZsig=TMath::Sqrt(radi);
268 goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
269 cout<<"Simmetria "<<TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))<<endl;
270 n1 = static_cast<Int_t>((fZFound-fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
271 n2 = static_cast<Int_t>((fZFound+fZCombc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
275 Warning("FindVertexForCurrentEvent","The procedure dows not converge\n");
278 cout<<"Numer of Iterations "<<niter<<endl<<endl;
279 fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
280 fCurrentVertex->SetTitle("vertexer: B");
282 return fCurrentVertex;
285 //_____________________________________________________________________
286 void AliITSVertexerZ::ResetHistograms(){
287 // delete TH1 data members
288 if(fZCombc)delete fZCombc;
289 if(fZCombf)delete fZCombf;
294 //______________________________________________________________________
295 void AliITSVertexerZ::FindVertices(){
296 // computes the vertices of the events in the range FirstEvent - LastEvent
297 AliRunLoader *rl = AliRunLoader::GetRunLoader();
298 AliITSLoader* itsLoader = (AliITSLoader*) rl->GetLoader("ITSLoader");
299 itsLoader->ReloadRecPoints();
300 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
301 cout<<"Processing event "<<i<<endl;
303 FindVertexForCurrentEvent(i);
305 WriteCurrentVertex();
309 cout<<"Vertex not found for event "<<i<<endl;
310 cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
316 //________________________________________________________
317 void AliITSVertexerZ::PrintStatus() const {
318 // Print current status
319 cout <<"=======================================================\n";
320 cout <<" First layer first and last modules: "<<fFirstL1<<", ";
321 cout <<fLastL1<<endl;
322 cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
323 cout <<fLastL2<<endl;
324 cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
325 cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
326 cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
327 cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
328 cout <<" Debug flag: "<<fDebug<<endl;
329 cout <<"First event to be processed "<<fFirstEvent;
330 cout <<"\n Last event to be processed "<<fLastEvent<<endl;
333 cout <<"=======================================================\n";