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