VertexerZ updated: no histograms for Z finding, new definition of errors
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerZ.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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<TBranch.h>
17 #include<TClonesArray.h>
18 #include<TFile.h>
19 #include<TH1.h>
20 #include <TString.h>
21 #include<TTree.h>
22 #include "AliITSLoader.h"
23 #include "AliITSgeom.h"
24 #include "AliITSDetTypeRec.h"
25 #include "AliITSRecPoint.h"
26 #include "AliITSZPoint.h"
27 #include "AliHeader.h"
28 #include "AliGenEventHeader.h"
29
30 /////////////////////////////////////////////////////////////////
31 // this class implements a fast method to determine
32 // the Z coordinate of the primary vertex
33 // for p-p collisions it seems to give comparable or better results
34 // with respect to what obtained with AliITSVertexerPPZ
35 // It can be used successfully with Pb-Pb collisions
36 ////////////////////////////////////////////////////////////////
37
38 ClassImp(AliITSVertexerZ)
39
40
41
42 //______________________________________________________________________
43 AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer(),
44 fFirstL1(0),
45 fLastL1(0),
46 fFirstL2(0),
47 fLastL2(0),
48 fDiffPhiMax(0),
49 fX0(0.),
50 fY0(0.),
51 fZFound(0),
52 fZsig(0.),
53 fZCombc(0),
54 fZCombf(0),
55 fLowLim(0.),
56 fHighLim(0.),
57 fStepCoarse(0),
58 fStepFine(0),
59 fTolerance(0.),
60 fMaxIter(0),
61 fWindowWidth(0) {
62   // Default constructor
63   SetDiffPhiMax();
64   SetFirstLayerModules();
65   SetSecondLayerModules();
66   SetLowLimit();
67   SetHighLimit();
68   SetBinWidthCoarse();
69   SetTolerance();
70   SetPPsetting();
71   ConfigIterations();
72   SetWindowWidth();
73 }
74
75 //______________________________________________________________________
76 AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn),
77 fFirstL1(0),
78 fLastL1(0),
79 fFirstL2(0),
80 fLastL2(0),
81 fDiffPhiMax(0),
82 fX0(x0),
83 fY0(y0),
84 fZFound(0),
85 fZsig(0.),
86 fZCombc(0),
87 fZCombf(0),
88 fLowLim(0.),
89 fHighLim(0.),
90 fStepCoarse(0),
91 fStepFine(0),
92 fTolerance(0.),
93 fMaxIter(0),
94 fWindowWidth(0) {
95   // Standard Constructor
96   SetDiffPhiMax();
97   SetFirstLayerModules();
98   SetSecondLayerModules();
99   SetBinWidthFine();
100   SetLowLimit();
101   SetHighLimit();
102   SetBinWidthCoarse();
103   SetTolerance();
104   SetPPsetting();
105   ConfigIterations();
106   SetWindowWidth();
107
108 }
109
110 //______________________________________________________________________
111 AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr),
112 fFirstL1(vtxr.fFirstL1),
113 fLastL1(vtxr.fLastL1),
114 fFirstL2(vtxr.fFirstL2),
115 fLastL2(vtxr.fLastL2),
116 fDiffPhiMax(vtxr.fDiffPhiMax),
117 fX0(vtxr.fX0),
118 fY0(vtxr.fY0),
119 fZFound(vtxr.fZFound),
120 fZsig(vtxr.fZsig),
121 fZCombc(vtxr.fZCombc),
122 fZCombf(vtxr.fZCombf),
123 fLowLim(vtxr.fLowLim),
124 fHighLim(vtxr.fHighLim),
125 fStepCoarse(vtxr.fStepCoarse),
126 fStepFine(vtxr.fStepFine),
127 fTolerance(vtxr.fTolerance),
128 fMaxIter(vtxr.fMaxIter),
129 fWindowWidth(vtxr.fWindowWidth){
130   // Copy constructor
131
132 }
133
134 //______________________________________________________________________
135 AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ&  vtxr ){
136   // Assignment operator
137
138   this->~AliITSVertexerZ();
139   new(this) AliITSVertexerZ(vtxr);
140   return *this;
141 }
142
143 //______________________________________________________________________
144 AliITSVertexerZ::~AliITSVertexerZ() {
145   // Destructor
146   delete fZCombc;
147 }
148
149 //______________________________________________________________________
150 void AliITSVertexerZ::ConfigIterations(Int_t noiter,Float_t *ptr){
151   // configure the iterative procedure to gain efficiency for
152   // pp events with very low multiplicity
153   Float_t defaults[5]={0.05,0.1,0.2,0.3,0.5};
154   fMaxIter=noiter;
155   if(noiter>5){
156     Error("ConfigIterations","Maximum number of iterations is 5\n");
157     fMaxIter=5;
158   }
159   for(Int_t j=0;j<5;j++)fPhiDiffIter[j]=defaults[j];
160   if(ptr)for(Int_t j=0;j<fMaxIter;j++)fPhiDiffIter[j]=ptr[j];
161 }
162
163 //______________________________________________________________________
164 Int_t AliITSVertexerZ::GetPeakRegion(TH1F*h, Int_t &binmin, Int_t &binmax) const {
165   // Finds a region around a peak in the Z histogram
166   // Case of 2 peaks is treated 
167   Int_t imax=h->GetNbinsX();
168   Float_t maxval=0;
169   Int_t bi1=h->GetMaximumBin();
170   Int_t bi2=0;
171   for(Int_t i=imax;i>=1;i--){
172     if(h->GetBinContent(i)>maxval){
173       maxval=h->GetBinContent(i);
174       bi2=i;
175     }
176   }
177   Int_t npeaks=0;
178
179   if(bi1==bi2){
180     binmin=bi1-3;
181     binmax=bi1+3;
182     npeaks=1;
183   }else{
184     TH1F *copy = new TH1F(*h);
185     copy->SetBinContent(bi1,0.);
186     copy->SetBinContent(bi2,0.);
187     Int_t l1=TMath::Max(bi1-3,1);
188     Int_t l2=TMath::Min(bi1+3,h->GetNbinsX());
189     Float_t cont1=copy->Integral(l1,l2);
190     Int_t ll1=TMath::Max(bi2-3,1);
191     Int_t ll2=TMath::Min(bi2+3,h->GetNbinsX());
192     Float_t cont2=copy->Integral(ll1,ll2);
193     if(cont1>cont2){
194       binmin=l1;
195       binmax=l2;
196       npeaks=1;
197     }
198     if(cont2>cont1){
199       binmin=ll1;
200       binmax=ll2;
201       npeaks=1;
202     }
203     if(cont1==cont2){
204       binmin=l1;
205       binmax=ll2;
206       if(bi2-bi1==1) npeaks=1;
207       else npeaks=2;
208     }  
209     delete copy;
210   }    
211   return npeaks;
212 }
213 //______________________________________________________________________
214 AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
215   // Defines the AliESDVertex for the current event
216   VertexZFinder(evnumber);
217   Int_t ntrackl=0;
218   for(Int_t iteraz=0;iteraz<fMaxIter;iteraz++){
219     if(fCurrentVertex) ntrackl=fCurrentVertex->GetNContributors();
220     if(!fCurrentVertex || ntrackl==0 || ntrackl==-1){
221       Float_t diffPhiMaxOrig=fDiffPhiMax;
222       fDiffPhiMax=GetPhiMaxIter(iteraz);
223       VertexZFinder(evnumber);
224       fDiffPhiMax=diffPhiMaxOrig;
225     }
226   }
227   FindMultiplicity(evnumber);
228   return fCurrentVertex;
229 }  
230
231
232
233
234 //______________________________________________________________________
235 void AliITSVertexerZ::VertexZFinder(Int_t evnumber){
236   // Defines the AliESDVertex for the current event
237   fCurrentVertex = 0;
238   AliRunLoader *rl =AliRunLoader::GetRunLoader();
239   AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
240   AliITSgeom* geom = itsLoader->GetITSgeom();
241   itsLoader->LoadRecPoints();
242   rl->GetEvent(evnumber);
243
244   AliITSDetTypeRec detTypeRec;
245
246   TTree *tR = itsLoader->TreeR();
247   detTypeRec.SetTreeAddressR(tR);
248   TClonesArray *itsRec  = 0;
249   // lc and gc are local and global coordinates for layer 1
250   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
251   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
252   // lc2 and gc2 are local and global coordinates for layer 2
253   Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
254   Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
255
256   itsRec = detTypeRec.RecPoints();
257   TBranch *branch;
258   branch = tR->GetBranch("ITSRecPoints");
259
260   Int_t nrpL1 = 0;
261   Int_t nrpL2 = 0;
262   // By default fFirstL1=0 and fLastL1=79
263   // This loop counts the number of recpoints on layer1 (central modules)
264   for(Int_t module= fFirstL1; module<=fLastL1;module++){
265     // Keep only central modules
266     //    if(module%4==0 || module%4==3)continue;
267     //   cout<<"Procesing module "<<module<<" ";
268     branch->GetEvent(module);
269     //    cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
270     nrpL1+= itsRec->GetEntries();
271     detTypeRec.ResetRecPoints();
272   }
273   //By default fFirstL2=80 and fLastL2=239
274   //This loop counts the number of RP on layer 2
275   for(Int_t module= fFirstL2; module<=fLastL2;module++){
276     branch->GetEvent(module);
277     nrpL2+= itsRec->GetEntries();
278     detTypeRec.ResetRecPoints();
279   }
280   if(nrpL1 == 0 || nrpL2 == 0){
281     ResetHistograms();
282     itsLoader->UnloadRecPoints();
283     fCurrentVertex = new AliESDVertex(0.,5.3,-2);
284     return;
285   }
286   // The vertex finding is attempted only if the number of RP is !=0 on
287   // both layers
288   Float_t *xc1 = new Float_t [nrpL1]; // coordinates of the L1 Recpoints
289   Float_t *yc1 = new Float_t [nrpL1];
290   Float_t *zc1 = new Float_t [nrpL1];
291   Float_t *phi1 = new Float_t [nrpL1];
292   Float_t *err1 = new Float_t [nrpL1];
293   Float_t *xc2 = new Float_t [nrpL2]; // coordinates of the L1 Recpoints
294   Float_t *yc2 = new Float_t [nrpL2];
295   Float_t *zc2 = new Float_t [nrpL2];
296   Float_t *phi2 = new Float_t [nrpL2];
297   Float_t *err2 = new Float_t [nrpL2];
298   Int_t ind = 0;// running index for RP
299   // Force a coarse bin size of 200 microns if the number of clusters on layer 2
300   // is low
301   if(nrpL2<fPPsetting[0])SetBinWidthCoarse(fPPsetting[1]);
302   // By default nbincoarse=(10+10)/0.01=2000
303   Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
304   if(fZCombc)delete fZCombc;
305   fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
306
307   // Loop on modules of layer 1 
308
309   for(Int_t module= fFirstL1; module<=fLastL1;module++){
310     //    if(module%4==0 || module%4==3)continue;
311     branch->GetEvent(module);
312     Int_t nrecp1 = itsRec->GetEntries();
313     for(Int_t j=0;j<nrecp1;j++){
314       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
315       // Local coordinates of this recpoint
316       lc[0]=recp->GetDetLocalX();
317       lc[2]=recp->GetDetLocalZ();
318       geom->LtoG(module,lc,gc);
319       // Global coordinates of this recpoints
320       gc[0]-=fX0; // Possible beam offset in the bending plane
321       gc[1]-=fY0; //   "               "
322       xc1[ind]=gc[0];
323       yc1[ind]=gc[1];
324       zc1[ind]=gc[2];
325       // azimuthal angle is computed in the interval 0 --> 2*pi
326       phi1[ind] = TMath::ATan2(gc[1],gc[0]);
327       if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
328       err1[ind]=recp->GetSigmaZ2();
329       ind++;
330     }
331     detTypeRec.ResetRecPoints();
332   }
333   ind = 0; // the running index is reset for Layer 2
334   for(Int_t module= fFirstL2; module<=fLastL2;module++){
335     branch->GetEvent(module);
336     Int_t nrecp2 = itsRec->GetEntries();
337     for(Int_t j=0;j<nrecp2;j++){
338       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
339       lc[0]=recp->GetDetLocalX();
340       lc[2]=recp->GetDetLocalZ();
341       geom->LtoG(module,lc,gc);
342       gc[0]-=fX0;
343       gc[1]-=fY0;
344       xc2[ind]=gc[0];
345       yc2[ind]=gc[1];
346       zc2[ind]=gc[2];
347       phi2[ind] = TMath::ATan2(gc[1],gc[0]);
348       if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
349       err2[ind]=recp->GetSigmaZ2();
350       ind++;
351     }
352     detTypeRec.ResetRecPoints();
353   }
354  
355 /* Test the ffect of mutiple scatternig on error. Negligible
356   // Multiple scattering
357   Float_t beta=1.,pmed=0.875; //pmed=875 MeV (for tracks with dphi<0.01 rad)
358   Float_t beta2=beta*beta;
359   Float_t p2=pmed*pmed;
360   Float_t rBP=3; //Beam Pipe radius = 3cm
361   Float_t dBP=0.08/35.3; // 800 um of Be
362   Float_t dL1=0.01; //approx. 1% of radiation length  
363   Float_t theta2BP=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dBP);
364   Float_t theta2L1=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(dL1);
365 */
366   TClonesArray *points = new TClonesArray("AliITSZPoint",nrpL1*nrpL2);
367   TClonesArray &pts = *points;
368   Int_t nopoints =0;
369   for(Int_t i=0;i<nrpL1;i++){ // loop on L1 RP
370     Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]); // radius L1 RP
371     for(Int_t j=0;j<nrpL2;j++){ // loop on L2 RP
372       Float_t diff = TMath::Abs(phi2[j]-phi1[i]); // diff in azimuth
373       if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; //diff<pi
374       if(diff<fDiffPhiMax){ // cut on 10 milliradians by def.
375         Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]); // radius L2 RP
376 //      Float_t tgth=(zc2[j]-zc1[i])/(r2-r1); // slope
377         Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1); //Z @ null radius
378         Float_t ezr0q=(r2*r2*err1[i]+r1*r1*err2[j])/(r2-r1)/(r2-r1); //error on Z @ null radius
379         /*
380         ezr0q+=r1*r1*(1+tgth*tgth)*theta2L1/2; // multiple scattering in layer 1
381         ezr0q+=rBP*rBP*(1+tgth*tgth)*theta2BP/2; // multiple scattering in beam pipe
382         */
383         new(pts[nopoints++])AliITSZPoint(zr0,ezr0q);
384
385         fZCombc->Fill(zr0);
386       }
387     }
388   }
389   delete [] xc1;
390   delete [] yc1;
391   delete [] zc1;
392   delete [] phi1;
393   delete [] err1;
394   delete [] xc2;
395   delete [] yc2;
396   delete [] zc2;
397   delete [] phi2;
398   delete [] err2;
399
400   points->Sort();
401
402   Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1);
403   if(contents<1.){
404     //    Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
405     ResetHistograms();
406     itsLoader->UnloadRecPoints();
407     fCurrentVertex = new AliESDVertex(0.,5.3,-1);
408     return;
409   }
410
411   TH1F *hc = fZCombc;
412
413   
414   if(hc->GetBinContent(hc->GetMaximumBin())<3)hc->Rebin(3);
415   Int_t binmin,binmax;
416   Int_t nPeaks=GetPeakRegion(hc,binmin,binmax);   
417   if(nPeaks==2)AliWarning("2 peaks found");
418   Float_t zm =0.;
419   Float_t ezm =0.;
420   Float_t lim1 = hc->GetBinLowEdge(binmin);
421   Float_t lim2 = hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax);
422
423   if(nPeaks ==1 && (lim2-lim1)<fWindowWidth){
424     Float_t c=(lim1+lim2)/2.;
425     lim1=c-fWindowWidth/2.;
426     lim2=c+fWindowWidth/2.;
427   }
428   Int_t niter = 0, ncontr=0;
429   do {
430     // symmetrization
431     if(zm  !=0.){
432       Float_t semilarg=TMath::Min((lim2-zm),(zm-lim1));
433       lim1=zm - semilarg;
434       lim2=zm + semilarg;
435     }
436
437     zm=0.;
438     ezm=0.;
439     ncontr=0;
440     for(Int_t i =0; i<points->GetEntriesFast(); i++){
441       AliITSZPoint* p=(AliITSZPoint*)points->UncheckedAt(i);
442       if(p->GetZ()>lim1 && p->GetZ()<lim2){
443         Float_t deno = p->GetErrZ();
444         zm+=p->GetZ()/deno;
445         ezm+=1./deno;
446         ncontr++;
447       }
448     }
449     zm/=ezm;
450     ezm=TMath::Sqrt(1./ezm);
451     niter++;
452   } while(niter<10 && TMath::Abs((zm-lim1)-(lim2-zm))>fTolerance);
453   fCurrentVertex = new AliESDVertex(zm,ezm,ncontr);
454   fCurrentVertex->SetTitle("vertexer: B");
455   ResetHistograms();
456   itsLoader->UnloadRecPoints();
457   return;
458 }
459
460 //_____________________________________________________________________
461 void AliITSVertexerZ::ResetHistograms(){
462   // delete TH1 data members
463   if(fZCombc)delete fZCombc;
464   fZCombc = 0;
465 }
466
467 //______________________________________________________________________
468 void AliITSVertexerZ::FindVertices(){
469   // computes the vertices of the events in the range FirstEvent - LastEvent
470   AliRunLoader *rl = AliRunLoader::GetRunLoader();
471   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
472   itsLoader->ReloadRecPoints();
473   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
474     //  cout<<"Processing event "<<i<<endl;
475     rl->GetEvent(i);
476     FindVertexForCurrentEvent(i);
477     if(fCurrentVertex){
478       WriteCurrentVertex();
479     }
480   }
481 }
482
483 //________________________________________________________
484 void AliITSVertexerZ::PrintStatus() const {
485   // Print current status
486   cout <<"=======================================================\n";
487   cout <<" First layer first and last modules: "<<fFirstL1<<", ";
488   cout <<fLastL1<<endl;
489   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
490   cout <<fLastL2<<endl;
491   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
492   cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
493   cout <<"Bin sizes for coarse z histos "<<fStepCoarse<<endl;
494   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
495   cout <<" Debug flag: "<<fDebug<<endl;
496   cout <<"First event to be processed "<<fFirstEvent;
497   cout <<"\n Last event to be processed "<<fLastEvent<<endl;
498   if(fZCombc){
499     cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl;
500   }
501   else{
502     cout<<"fZCombc does not exist\n";
503   }
504  
505   cout <<"=======================================================\n";
506 }
507