Print removed
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerPPZ.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 <AliITSVertexerPPZ.h>
16 #include <Riostream.h>
17 #include <TArrayF.h>
18 #include <TDirectory.h>
19 #include <TH1F.h>
20 #include <TMath.h>
21 #include <TObjArray.h>
22 #include <TTree.h>
23 #include <TClonesArray.h>
24 #include "AliRun.h"
25 #include "AliITS.h"
26 #include "AliITSgeom.h"
27 #include "AliITSLoader.h"
28 #include "AliITSRecPoint.h"
29 /////////////////////////////////////////////////////////////////////////
30 //                                                                     //
31 // This class is intended to compute the Z coordinate of the           //
32 // primary vertex for p-p interactions, or in general when the         //
33 // number of secondaries is too slow to use the class                  //
34 // AliITSVertexerIons, which in turn is optimized for A-A collsions    //
35 // Origin: masera@to.infn.it      9/12/2002                            //
36 //                                                                     //
37 /////////////////////////////////////////////////////////////////////////
38 ClassImp(AliITSVertexerPPZ)
39
40
41
42 //______________________________________________________________________
43 AliITSVertexerPPZ::AliITSVertexerPPZ():AliITSVertexer() {
44   // Default Constructor
45
46   SetDiffPhiMax(0);
47   fX0 = 0.;
48   fY0 = 0.;
49   SetFirstLayerModules(0);
50   SetSecondLayerModules(0);
51   fZFound = 0;
52   fZsig = 0.;
53   fITS = 0;
54   SetWindow(0);
55 }
56
57 AliITSVertexerPPZ::AliITSVertexerPPZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
58   // Standard constructor
59   SetDiffPhiMax();
60   fX0 = x0;
61   fY0 = y0;
62   SetFirstLayerModules();
63   SetSecondLayerModules();
64   fZFound = 0;
65   fZsig = 0.;
66   fITS = 0;
67   SetWindow();
68 }
69
70 //______________________________________________________________________
71 AliITSVertexerPPZ::~AliITSVertexerPPZ() {
72   // Default Destructor
73   fITS = 0;
74 }
75
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
79   fZFound=0;
80   fZsig=0;
81   Int_t N=0;
82   Int_t NbinNotZero=0;
83   Float_t totst = 0.;
84   Float_t totst2 = 0.;
85   Float_t curz = 0.;
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);
89     totst+=cont;
90     totst2+=cont*cont;
91     N++;
92     if(cont!=0)NbinNotZero++;
93   }
94   if(NbinNotZero==0){fZFound=-100; fZsig=-100; return;}
95   if(NbinNotZero==1){
96     fZFound = curz;
97     fZsig=0;
98     fCurrentVertex = new AliESDVertex(fZFound,fZsig,NbinNotZero);
99     return;
100   }
101   Float_t errsq = totst2/(N-1)-totst*totst/N/(N-1);
102   if(errsq>=0){
103   totst2=TMath::Sqrt(errsq);
104   }
105   else {
106     Error("EvalZ","Negative variance: %d - %12.7f - %12.7f",N,totst2,totst);
107     fZFound=-100; 
108     fZsig=-100;
109     return;
110   }
111   totst /= N;
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);
116   Float_t valm = 0.;
117   Float_t zmax = 0.;
118   for(Int_t i=1;i<=sepa;i++){
119     Float_t cont=hist->GetBinContent(i);
120     if(cont>valm){
121       valm = cont;
122       zmax = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(1);
123     }
124     if(cont>cut){
125       curz=hist->GetBinLowEdge(i);
126       if(curz<val1)val1=curz;
127       if(curz>val2)val2=curz;
128     }
129   }
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";
135   }
136   if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<" "<<val2-val1<<endl;
137   fZFound=0;
138   fZsig=0;
139   N=0;
140   for(Int_t i=0; i<ncoinc; i++){
141     Float_t z=zval->At(i);
142     if(z<val1)continue;
143     if(z>val2)continue;
144    
145     fZFound+=z;
146     fZsig+=z*z;
147     /*   weights defined by the curvature
148     Float_t wei = 1./curv->At(i);
149     fZFound+=z*wei;
150     fZsig+=wei;
151     */
152     N++;
153   }
154   if(N<1){fZFound=-110; fZsig=-110; return;}
155   if(N==1){
156     fZsig = 0;
157     fCurrentVertex = new AliESDVertex(fZFound,fZsig,N);
158     return;
159   }
160   errsq = (fZsig/(N-1)-fZFound*fZFound/N/(N-1))/N;
161   if(errsq>=0.){
162     fZsig=TMath::Sqrt(errsq);
163   }
164   else {
165     Error("evalZ","Negative variance: %d - %12.7f %12.7f",N,fZsig,fZFound);
166     fZsig=0;
167   }
168   fZFound=fZFound/N;
169   /* weights defined by the curvature
170   fZsig=1./fZsig;
171   fZFound*=fZsig;
172   fZsig = TMath::Sqrt(fZsig);
173   */
174   fCurrentVertex = new AliESDVertex(fZFound,fZsig,N);
175 }
176
177 //______________________________________________________________________
178 AliESDVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
179   // Defines the AliESDVertex for the current event
180   fCurrentVertex = 0;
181   fZFound = -999;
182   fZsig = -999;
183   AliRunLoader *rl =AliRunLoader::GetRunLoader();
184   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
185   if(!fITS)  {
186     fITS = (AliITS*)gAlice->GetModule("ITS");
187     if(!fITS) {
188       Error("FindVertexForCurrentEvent","AliITS object was not found");
189       return fCurrentVertex;
190     }
191   }
192   fITS->SetTreeAddress();
193   AliITSgeom *geom = fITS->GetITSgeom();
194   if(!geom) {
195     Error("FindVertexForCurrentEvent","ITS geometry is not defined");
196     return fCurrentVertex;
197   }
198   TTree *tR=0;
199   TClonesArray *itsRec  = 0;
200   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
201   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
202   Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
203   Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
204
205   tR = ITSloader->TreeR();
206   if(!tR){
207     Error("FindVertexForCurrentEvent","TreeR not found");
208     return fCurrentVertex;
209   }
210   itsRec  = fITS->RecPoints();  // uses slow points or fast points if slow are
211   // missing
212   TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
213   TBranch *branch;
214   if(fUseV2Clusters){
215     branch = tR->GetBranch("Clusters");
216     branch->SetAddress(&clusters);
217   }
218   else {
219     branch = tR->GetBranch("ITSRecPoints");
220     if(!branch){ 
221       branch = tR->GetBranch("ITSRecPointsF");
222     }
223   }
224   if(!branch){
225    Error("FindVertexForCurrentEvent","branch for ITS rec points not found");
226    return fCurrentVertex;
227   }
228   Float_t zave=0;
229   Float_t rmszav=0;
230   Float_t zave2=0;
231   Int_t firipixe=0;
232   for(Int_t module= fFirstL1; module<=fLastL1;module++){
233     branch->GetEvent(module);
234     if(fUseV2Clusters){
235       Clusters2RecPoints(clusters,module,itsRec);
236     }
237     Int_t nrecp1 = itsRec->GetEntries();
238     for(Int_t i=0; i<nrecp1;i++){
239       AliITSRecPoint *current = (AliITSRecPoint*)itsRec->At(i);
240       lc[0]=current->GetX();
241       lc[2]=current->GetZ();
242       geom->LtoG(module,lc,gc);
243       zave+=gc[2];
244       zave2+=gc[2]*gc[2];
245       firipixe++;
246     }
247     fITS->ResetRecPoints();
248   }
249   if(firipixe>1){
250     rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
251     zave=zave/firipixe;
252     if(fDebug>1)cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++\n Number of firing pixels: "<<firipixe<<endl;
253   }
254   else {
255     fZFound = -200;
256     fZsig = -200;
257     Warning("FindVertexForCurrentEvent","No rec points on first layer for this event");
258     return fCurrentVertex;
259   }
260   Float_t zlim1=zave-rmszav;
261   Float_t zlim2=zave+rmszav;
262   Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
263   zlim2=zlim1 + sepa/10.;
264   TH1F *zvdis = new TH1F("z_ev","zv distr",sepa,zlim1,zlim2);
265   if(fDebug>0){
266     cout<<"Z limits: "<<zlim1<<" "<<zlim2<<"; Bins= "<<sepa<<endl;
267     cout<<"Bin width: "<<zvdis->GetBinWidth(1)*10000.<<" micron\n";
268   }
269   Int_t sizarr=100;
270   TArrayF *zval = new TArrayF(sizarr);
271   //  TArrayF *curv = new TArrayF(sizarr);
272   Int_t ncoinc=0;
273   for(Int_t module= fFirstL1; module<=fLastL1;module++){
274     if(fDebug>0)cout<<"processing module   "<<module<<"                  \r";
275     branch->GetEvent(module);
276     if(fUseV2Clusters){
277       Clusters2RecPoints(clusters,module,itsRec);
278     }
279     Int_t nrecp1 = itsRec->GetEntries();
280     TObjArray *poiL1 = new TObjArray(nrecp1);
281     for(Int_t i=0; i<nrecp1;i++)poiL1->AddAt(itsRec->At(i),i);
282     fITS->ResetRecPoints();
283     for(Int_t i=0; i<nrecp1;i++){
284       AliITSRecPoint *current = (AliITSRecPoint*)poiL1->At(i);
285       lc[0]=current->GetX();
286       lc[2]=current->GetZ();
287       geom->LtoG(module,lc,gc);
288       gc[0]-=fX0;
289       gc[1]-=fY0;
290       Float_t r1=TMath::Sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
291       Float_t phi1 = TMath::ATan2(gc[1],gc[0]);
292       if(phi1<0)phi1=2*TMath::Pi()+phi1;
293       if(fDebug>1)cout<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1<<"     \n";
294       for(Int_t modul2=fFirstL2; modul2<=fLastL2; modul2++){
295         branch->GetEvent(modul2);
296         if(fUseV2Clusters){
297           Clusters2RecPoints(clusters,modul2,itsRec);
298         }
299         Int_t nrecp2 = itsRec->GetEntries();
300         for(Int_t j=0; j<nrecp2;j++){
301           AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
302           lc2[0]=recp->GetX();
303           lc2[2]=recp->GetZ();
304           geom->LtoG(modul2,lc2,gc2);
305           gc2[0]-=fX0;
306           gc2[1]-=fY0;
307           Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
308           Float_t zr0=(r2*gc[2]-r1*gc2[2])/(r2-r1);
309           Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
310           if(phi2<0)phi2=2.*TMath::Pi()+phi2;
311           Float_t diff = TMath::Abs(phi2-phi1);
312           if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
313           if(zr0>zlim1 && zr0<zlim2){
314             if(diff<fDiffPhiMax ){
315               zvdis->Fill(zr0);
316               zval->AddAt(zr0,ncoinc);
317               /* uncomment these lines to use curvature as a weight
318               Float_t cu = Curv(0.,0.,gc[0],gc[1],gc2[0],gc[1]);
319               curv->AddAt(cu,ncoinc);
320               */
321               ncoinc++;
322               if(ncoinc==(sizarr-1)){
323                 sizarr+=100;
324                 zval->Set(sizarr);
325                 //uncomment next line to use curvature as weight
326                 //              curv->Set(sizarr);
327               }
328             }
329           }
330         }
331         fITS->ResetRecPoints();
332       }
333     }
334     delete poiL1;
335   }         // loop on modules
336   if(fDebug>0){
337     cout<<endl<<"Number of coincidences = "<<ncoinc<<endl;
338   }
339   //  EvalZ(zvdis,sepa,ncoinc,zval,curv);
340   EvalZ(zvdis,sepa,ncoinc,zval);
341   delete zvdis;
342   delete zval;
343   //  delete curv;
344   if(fCurrentVertex){
345     char name[30];
346     sprintf(name,"Vertex_%d",evnumber);
347     //    fCurrentVertex->SetName(name);
348     fCurrentVertex->SetTitle("vertexer: PPZ");
349   }
350   return fCurrentVertex;
351 }
352
353 //______________________________________________________________________
354 void AliITSVertexerPPZ::FindVertices(){
355   // computes the vertices of the events in the range FirstEvent - LastEvent
356   AliRunLoader *rl = AliRunLoader::GetRunLoader();
357   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
358   ITSloader->ReloadRecPoints();
359   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
360     cout<<"Processing event "<<i<<endl;
361     rl->GetEvent(i);
362     FindVertexForCurrentEvent(i);
363     if(fCurrentVertex){
364       WriteCurrentVertex();
365     }
366     else {
367       if(fDebug>0){
368         cout<<"Vertex not found for event "<<i<<endl;
369         cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
370       }
371     }
372   }
373 }
374
375 //________________________________________________________
376 void AliITSVertexerPPZ::PrintStatus() const {
377   // Print current status
378   cout <<"=======================================================\n";
379   cout <<" First layer first and last modules: "<<fFirstL1<<", ";
380   cout <<fLastL1<<endl;
381   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
382   cout <<fLastL2<<endl;
383   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
384   cout <<" Window for Z search: "<<fWindow<<endl;
385   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
386   cout <<" Debug flag: "<<fDebug<<endl;
387   cout<<"First event to be processed "<<fFirstEvent;
388   cout<<"\n Last event to be processed "<<fLastEvent<<endl;
389 }
390
391 //________________________________________________________
392 Float_t AliITSVertexerPPZ::Curv(Double_t x1,Double_t y1,
393                    Double_t x2,Double_t y2,
394                    Double_t x3,Double_t y3) 
395 {
396   //-----------------------------------------------------------------
397   // Initial approximation of the track curvature (Y. Belikov) squared
398   //-----------------------------------------------------------------
399   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
400   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
401                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
402   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
403                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
404
405   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
406
407   Float_t val = (Float_t)(( -xr*yr/TMath::Sqrt(xr*xr+yr*yr))
408                           *( -xr*yr/TMath::Sqrt(xr*xr+yr*yr)));
409   return val; 
410 }