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 <AliITSVertexerPPZ.h>
16 #include <Riostream.h>
18 #include <TDirectory.h>
21 #include <TObjArray.h>
23 #include <TClonesArray.h>
26 #include "AliITSgeom.h"
27 #include "AliITSRecPoint.h"
28 /////////////////////////////////////////////////////////////////////////
30 // This class is intended to compute the Z coordinate of the //
31 // primary vertex for p-p interactions, or in general when the //
32 // number of secondaries is too slow to use the class //
33 // AliITSVertexerIons, which in turn is optimized for A-A collsions //
34 // Origin: masera@to.infn.it 9/12/2002 //
36 /////////////////////////////////////////////////////////////////////////
37 ClassImp(AliITSVertexerPPZ)
41 //______________________________________________________________________
42 AliITSVertexerPPZ::AliITSVertexerPPZ():AliITSVertexer() {
43 // Default Constructor
48 SetFirstLayerModules(0);
49 SetSecondLayerModules(0);
56 AliITSVertexerPPZ::AliITSVertexerPPZ(TFile *infile, TFile *outfile, Float_t x0, Float_t y0):AliITSVertexer(infile,outfile) {
57 // Standard constructor
58 if(!fInFile)Fatal("AliITSVertexerPPZ","No inputfile has been defined!");
62 SetFirstLayerModules();
63 SetSecondLayerModules();
70 //______________________________________________________________________
71 AliITSVertexerPPZ::~AliITSVertexerPPZ() {
76 //________________________________________________________
77 void AliITSVertexerPPZ::EvalZ(TH1F *hist,Int_t sepa, Int_t ncoinc, TArrayF *zval) {
78 Float_t DeltaVal = hist->GetBinWidth(1)*fWindow; // max window in Z for searching
86 for(Int_t i=1;i<=sepa;i++){
87 Float_t cont=hist->GetBinContent(i);
88 if(cont!=0)curz = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(i);
92 if(cont!=0)NbinNotZero++;
94 if(NbinNotZero==0){fZFound=-100; fZsig=-100; return;}
98 fCurrentVertex = new AliITSVertex(fZFound,fZsig,NbinNotZero);
101 Float_t errsq = totst2/(N-1)-totst*totst/N/(N-1);
103 totst2=TMath::Sqrt(errsq);
106 Error("EvalZ","Negative variance: %d - %12.7f - %12.7f",N,totst2,totst);
112 Float_t cut = totst+totst2*2.;
113 if(fDebug>1)cout<<"totst, totst2, cut: "<<totst<<", "<<totst2<<", "<<cut<<endl;
114 Float_t val1=hist->GetBinLowEdge(sepa);
115 Float_t val2=hist->GetBinLowEdge(1);
118 for(Int_t i=1;i<=sepa;i++){
119 Float_t cont=hist->GetBinContent(i);
122 zmax = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(1);
125 curz=hist->GetBinLowEdge(i);
126 if(curz<val1)val1=curz;
127 if(curz>val2)val2=curz;
130 val2+=hist->GetBinWidth(1);
131 if((val2-val1)>DeltaVal){
132 val1 = zmax-DeltaVal/2.;
133 val2 = zmax+DeltaVal/2.;
134 if(fDebug>0)cout<<"val1 and val2 recomputed\n";
136 if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<" "<<val2-val1<<endl;
140 for(Int_t i=0; i<ncoinc; i++){
141 Float_t z=zval->At(i);
147 /* weights defined by the curvature
148 Float_t wei = 1./curv->At(i);
154 if(N<1){fZFound=-110; fZsig=-110; return;}
157 fCurrentVertex = new AliITSVertex(fZFound,fZsig,N);
160 errsq = (fZsig/(N-1)-fZFound*fZFound/N/(N-1))/N;
162 fZsig=TMath::Sqrt(errsq);
165 Error("evalZ","Negative variance: %d - %12.7f %12.7f",N,fZsig,fZFound);
169 /* weights defined by the curvature
172 fZsig = TMath::Sqrt(fZsig);
174 fCurrentVertex = new AliITSVertex(fZFound,fZsig,N);
177 //______________________________________________________________________
178 AliITSVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
179 // Defines the AliITSVertex for the current event
183 if(!gAlice)gAlice = (AliRun*)fInFile->Get("gAlice");
185 Error("FindVertexForCurrentEvent","The AliRun object is not available - nothing done");
186 return fCurrentVertex;
189 fITS = (AliITS*)gAlice->GetModule("ITS");
191 Error("FindVertexForCurrentEvent","AliITS object was not found");
192 return fCurrentVertex;
195 AliITSgeom *geom = fITS->GetITSgeom();
197 Error("FindVertexForCurrentEvent","ITS geometry is not defined");
198 return fCurrentVertex;
201 TClonesArray *ITSrec = 0;
202 Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
203 Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
204 Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
205 Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
207 TR = gAlice->TreeR();
209 Error("FindVertexForCurrentEvent","TreeR not found");
210 return fCurrentVertex;
212 ITSrec = fITS->RecPoints(); // uses slow points or fast points if slow are
214 TBranch *branch = TR->GetBranch("ITSRecPoints");
216 branch = TR->GetBranch("ITSRecPointsF");
217 // Warning("FindVertexForCurrentEvent","Using Fast points");
220 Error("FindVertexForCurrentEvent","branch for ITS rec points not found");
221 return fCurrentVertex;
227 for(Int_t module= fFirstL1; module<=fLastL1;module++){
228 branch->GetEvent(module);
229 Int_t nrecp1 = ITSrec->GetEntries();
230 for(Int_t i=0; i<nrecp1;i++){
231 AliITSRecPoint *current = (AliITSRecPoint*)ITSrec->At(i);
232 lc[0]=current->GetX();
233 lc[2]=current->GetZ();
234 geom->LtoG(module,lc,gc);
239 fITS->ResetRecPoints();
242 rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
244 if(fDebug>1)cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++\n Number of firing pixels: "<<firipixe<<endl;
249 Warning("FindVertexForCurrentEvent","No rec points on first layer for this event");
250 return fCurrentVertex;
252 Float_t zlim1=zave-rmszav;
253 Float_t zlim2=zave+rmszav;
254 Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
255 zlim2=zlim1 + sepa/10.;
256 TH1F *zvdis = new TH1F("z_ev","zv distr",sepa,zlim1,zlim2);
258 cout<<"Z limits: "<<zlim1<<" "<<zlim2<<"; Bins= "<<sepa<<endl;
259 cout<<"Bin width: "<<zvdis->GetBinWidth(1)*10000.<<" micron\n";
262 TArrayF *zval = new TArrayF(sizarr);
263 // TArrayF *curv = new TArrayF(sizarr);
265 for(Int_t module= fFirstL1; module<=fLastL1;module++){
266 if(fDebug>0)cout<<"processing module "<<module<<" \r";
267 branch->GetEvent(module);
268 Int_t nrecp1 = ITSrec->GetEntries();
269 TObjArray *poiL1 = new TObjArray(nrecp1);
270 for(Int_t i=0; i<nrecp1;i++)poiL1->AddAt(ITSrec->At(i),i);
271 fITS->ResetRecPoints();
272 for(Int_t i=0; i<nrecp1;i++){
273 AliITSRecPoint *current = (AliITSRecPoint*)poiL1->At(i);
274 lc[0]=current->GetX();
275 lc[2]=current->GetZ();
276 geom->LtoG(module,lc,gc);
279 Float_t r1=TMath::Sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
280 Float_t phi1 = TMath::ATan2(gc[1],gc[0]);
281 if(phi1<0)phi1=2*TMath::Pi()+phi1;
282 if(fDebug>1)cout<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1<<" \n";
283 for(Int_t modul2=fFirstL2; modul2<=fLastL2; modul2++){
284 branch->GetEvent(modul2);
285 Int_t nrecp2 = ITSrec->GetEntries();
286 for(Int_t j=0; j<nrecp2;j++){
287 AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(j);
290 geom->LtoG(modul2,lc2,gc2);
293 Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
294 Float_t zr0=(r2*gc[2]-r1*gc2[2])/(r2-r1);
295 Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
296 if(phi2<0)phi2=2.*TMath::Pi()+phi2;
297 Float_t diff = TMath::Abs(phi2-phi1);
298 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
299 if(zr0>zlim1 && zr0<zlim2){
300 if(diff<fDiffPhiMax ){
302 zval->AddAt(zr0,ncoinc);
303 /* uncomment these lines to use curvature as a weight
304 Float_t cu = Curv(0.,0.,gc[0],gc[1],gc2[0],gc[1]);
305 curv->AddAt(cu,ncoinc);
308 if(ncoinc==(sizarr-1)){
311 //uncomment next line to use curvature as weight
312 // curv->Set(sizarr);
317 fITS->ResetRecPoints();
323 cout<<endl<<"Number of coincidences = "<<ncoinc<<endl;
325 // EvalZ(zvdis,sepa,ncoinc,zval,curv);
326 EvalZ(zvdis,sepa,ncoinc,zval);
332 sprintf(name,"Vertex_%d",evnumber);
333 fCurrentVertex->SetName(name);
334 fCurrentVertex->SetTitle("vertexer: PPZ");
336 return fCurrentVertex;
339 //______________________________________________________________________
340 void AliITSVertexerPPZ::FindVertices(){
341 // computes the vertices of the events in the range FirstEvent - LastEvent
343 Error("FindVertices","The output file is not available - nothing done");
346 if(!gAlice)gAlice = (AliRun*)fInFile->Get("gAlice");
348 Error("FindVertices","The AliRun object is not available - nothing done");
351 TDirectory *curdir = gDirectory;
353 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
355 FindVertexForCurrentEvent(i);
357 WriteCurrentVertex();
361 cout<<"Vertex not found for event "<<i<<endl;
362 cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
369 //________________________________________________________
370 void AliITSVertexerPPZ::PrintStatus() const {
371 // Print current status
372 cout <<"=======================================================\n";
373 cout <<" First layer first and last modules: "<<fFirstL1<<", ";
374 cout <<fLastL1<<endl;
375 cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
376 cout <<fLastL2<<endl;
377 cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
378 cout <<" Window for Z search: "<<fWindow<<endl;
379 cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
380 cout <<" Debug flag: "<<fDebug<<endl;
381 if(fInFile)cout <<" The input file is open\n";
382 if(!fInFile)cout <<"The input file is not open\n";
383 if(fOutFile)cout <<"The output file is open\n";
384 if(!fOutFile)cout <<"The output file is not open\n";
385 cout<<"First event to be processed "<<fFirstEvent;
386 cout<<"\n Last event to be processed "<<fLastEvent<<endl;
389 //________________________________________________________
390 Float_t AliITSVertexerPPZ::Curv(Double_t x1,Double_t y1,
391 Double_t x2,Double_t y2,
392 Double_t x3,Double_t y3)
394 //-----------------------------------------------------------------
395 // Initial approximation of the track curvature (Y. Belikov) squared
396 //-----------------------------------------------------------------
397 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
398 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
399 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
400 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
401 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
403 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
405 Float_t val = (Float_t)(( -xr*yr/TMath::Sqrt(xr*xr+yr*yr))
406 *( -xr*yr/TMath::Sqrt(xr*xr+yr*yr)));