Changes to obey coding conventions
[u/mrichter/AliRoot.git] / ITS / AliITSOnlineSDDInjectors.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-2009, 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 "AliITSOnlineSDDInjectors.h"
16 #include <TH2F.h>
17 #include <TGraphErrors.h>
18 #include <TMath.h>
19
20 ///////////////////////////////////////////////////////////////////
21 //                                                               //
22 // Implementation of the class used for SDD injector analysis    //
23 // Origin: F.Prino, Torino, prino@to.infn.it                     //
24 //                                                               //
25 ///////////////////////////////////////////////////////////////////
26
27 ClassImp(AliITSOnlineSDDInjectors)
28
29 const Float_t AliITSOnlineSDDInjectors::fgkSaturation = 1008.;
30 const Float_t AliITSOnlineSDDInjectors::fgkJitterTB = 8.;
31
32 //______________________________________________________________________
33   AliITSOnlineSDDInjectors::AliITSOnlineSDDInjectors():AliITSOnlineSDD(),fHisto(),fTbZero(0.),fParam(),fPolOrder(0),fMinDriftVel(0.),fMaxDriftVel(0.),fThreshold(0.)
34 {
35   // default constructor
36   SetMinDriftVel();
37   SetMaxDriftVel();
38   SetRangeLine1();
39   SetRangeLine2();
40   SetRangeLine3();
41   SetPositions();
42   SetPolOrder();
43   SetThreshold();
44 }
45 //______________________________________________________________________
46 AliITSOnlineSDDInjectors::AliITSOnlineSDDInjectors(Int_t mod, Int_t sid):AliITSOnlineSDD(mod,sid),fHisto(),fTbZero(0.),fParam(),fPolOrder(0),fMinDriftVel(0.),fMaxDriftVel(0.),fThreshold(0.)
47
48 // standard constructor
49   SetMinDriftVel();
50   SetMaxDriftVel();
51   SetRangeLine1();
52   SetRangeLine2();
53   SetRangeLine3();
54   SetPositions();
55   SetPolOrder();
56   SetThreshold();
57 }
58 //______________________________________________________________________
59 AliITSOnlineSDDInjectors::~AliITSOnlineSDDInjectors(){
60   // Destructor
61   if(fHisto) delete fHisto;  
62   if(fParam) delete [] fParam;
63 }
64 //______________________________________________________________________
65 void AliITSOnlineSDDInjectors::SetPositions(){
66   // 
67   Float_t kLinFromCenterUm[3]={31860.,17460.,660.};
68   Float_t kAnodeFromCenterUm=35085;
69   for(Int_t i=0;i<3;i++){
70     fPosition[i]=kAnodeFromCenterUm-kLinFromCenterUm[i];
71     fPosition[i]/=10000.; // from microns to cm
72   }
73 }
74 //______________________________________________________________________
75 void AliITSOnlineSDDInjectors::Reset(){
76   //
77   for(Int_t i=0;i<kNInjectors;i++){ 
78     fDriftVel[i]=0.;
79     fSigmaDriftVel[i]=0.;
80   }
81   for(Int_t i=0;i<kNInjectors;i++){
82     for(Int_t j=0;j<3;j++){
83       fGoodInj[i][j]=0;
84       fCentroid[i][j]=0.;
85       fRMSCentroid[i][j]=0.;
86     }
87   }
88 }
89 //______________________________________________________________________
90 void AliITSOnlineSDDInjectors::AnalyzeEvent(TH2F* his){
91   //
92   Reset();
93   fHisto=his;
94   FindGoodInjectors();
95   FindCentroids();
96   CalcTimeBinZero();
97   for(Int_t j=0;j<kNInjectors;j++) CalcDriftVelocity(j);
98   FitDriftVelocityVsAnode();
99 }
100 //______________________________________________________________________
101 TGraphErrors* AliITSOnlineSDDInjectors::GetLineGraph(Int_t jlin){
102   // 
103   Float_t x[4],y[4],ex[4],ey[4];
104   x[0]=0.;
105   ex[0]=0.;
106   y[0]=fTbZero;
107   ey[0]=0.;
108   for(Int_t i=0;i<3;i++){
109     x[i+1]=fPosition[i];
110     ex[i+1]=0.;
111     y[i+1]=fCentroid[jlin][i];
112     ey[i+1]=fRMSCentroid[jlin][i];
113   }
114   TGraphErrors *g=new TGraphErrors(4,x,y,ex,ey);
115   return g;
116 }
117 //______________________________________________________________________
118 Float_t AliITSOnlineSDDInjectors::GetDriftCoordinate(Float_t cAnode, Float_t cTimeBin){
119   //
120   Float_t vel=0;
121   for(Int_t i=0;i<=fPolOrder;i++) vel+=fParam[i]*TMath::Power(cAnode,(Float_t)i);
122   return vel*(cTimeBin-(fTbZero-fgkJitterTB))*25/1000.; 
123 }
124 //______________________________________________________________________
125 TGraphErrors* AliITSOnlineSDDInjectors::GetDriftVelocityGraph() const{
126   // 
127   Int_t ipt=0;
128   TGraphErrors *g=new TGraphErrors(0);
129   for(Int_t i=0;i<kNInjectors;i++){
130     if(fDriftVel[i]>0){ 
131       g->SetPoint(ipt,GetAnodeNumber(i),fDriftVel[i]);
132       g->SetPointError(ipt,0,fSigmaDriftVel[i]);
133       ipt++;
134     }
135   }
136   return g;
137 }
138 //______________________________________________________________________
139 void AliITSOnlineSDDInjectors::CalcTimeBinZero(){
140   //
141   Float_t tzero=0.,intCont=0.;
142   for(Int_t ian=0;ian<fgkNAnodes;ian++){
143     for(Int_t itb=1;itb<fTbMin[0];itb++){
144       Float_t cont=fHisto->GetBinContent(itb,ian+1);
145       if(cont>fThreshold){
146         tzero+=cont*float(itb);
147         intCont+=cont;
148       }
149     }
150   }
151   if(intCont>0) fTbZero=tzero/intCont;
152 }
153 //______________________________________________________________________
154 void AliITSOnlineSDDInjectors::FitDriftVelocityVsAnode(){
155   const Int_t nn=fPolOrder+1;
156   Float_t **mat = new Float_t*[nn];
157   for(Int_t i=0; i < nn; i++) mat[i] = new Float_t[nn];
158   Float_t *vect = new Float_t[nn];
159   for(Int_t k1=0;k1<nn;k1++){
160     vect[k1]=0;
161     for(Int_t k2=0;k2<nn;k2++){
162       mat[k1][k2]=0;
163       for(Int_t n=0; n<kNInjectors;n++){
164         Float_t x=(Float_t)GetAnodeNumber(n);
165         if(fDriftVel[n]>0) mat[k1][k2]+=TMath::Power(x,k1+k2)/TMath::Power(fSigmaDriftVel[n],2);
166       }
167     }
168   }
169   for(Int_t k1=0;k1<nn;k1++){
170     for(Int_t n=0; n<kNInjectors;n++){
171       Float_t x=(Float_t)GetAnodeNumber(n);
172       if(fDriftVel[n]>0) vect[k1]+=fDriftVel[n]*TMath::Power(x,k1)/TMath::Power(fSigmaDriftVel[n],2);
173     }
174   }
175   Int_t *iPivot = new Int_t[nn];
176   Int_t *indxR = new Int_t[nn];
177   Int_t *indxC = new Int_t[nn];
178   for(Int_t i=0;i<nn;i++) iPivot[i]=0;
179   Int_t iCol=-1,iRow=-1;
180   for(Int_t i=0;i<nn;i++){
181     Float_t big=0.;
182     for(Int_t j=0;j<nn;j++){
183       if(iPivot[j]!=1){
184         for(Int_t k=0;k<nn;k++){
185            if(iPivot[k]==0){
186              if(TMath::Abs(mat[j][k])>=big){
187                big=TMath::Abs(mat[j][k]);
188                iRow=j;
189                iCol=k;
190              }
191           }
192         }
193       }
194     }
195     iPivot[iCol]++;
196     Float_t aux;
197     if(iRow!=iCol){
198       for(Int_t l=0;l<nn;l++){
199         aux=mat[iRow][l];
200         mat[iRow][l]=mat[iCol][l];
201         mat[iCol][l]=aux;
202       }
203       aux=vect[iRow];
204       vect[iRow]=vect[iCol];
205       vect[iCol]=aux;
206     }
207     indxR[i]=iRow;
208     indxC[i]=iCol;
209     if(mat[iCol][iCol]==0) break;
210     Float_t pivinv=1./mat[iCol][iCol];
211     mat[iCol][iCol]=1;
212     for(Int_t l=0;l<nn;l++) mat[iCol][l]*=pivinv;
213     vect[iCol]*=pivinv;
214     for(Int_t m=0;m<nn;m++){
215       if(m!=iCol){
216         aux=mat[m][iCol];
217         mat[m][iCol]=0;
218         for(Int_t n=0;n<nn;n++) mat[m][n]-=mat[iCol][n]*aux;
219         vect[m]-=vect[iCol]*aux;
220       }
221     }    
222   }
223   delete [] iPivot;
224   delete [] indxR;
225   delete [] indxC;
226
227   if(fParam) delete [] fParam;
228   fParam=new Float_t[nn];
229   for(Int_t i=0; i<nn;i++)fParam[i]=vect[i];
230
231   for(Int_t i=0; i < nn; i++) delete [] mat[i];
232   delete [] mat;
233   delete [] vect;
234 }
235 //______________________________________________________________________
236 void AliITSOnlineSDDInjectors::CalcDriftVelocity(Int_t jlin){
237   // 
238   Float_t sumY=0,sumX=0,sumXX=0,sumYY=0.,sumXY=0,sumWEI=0.;
239   Int_t npt=0;
240   Float_t y[3],ey[3];
241   Float_t tzero=0,erry=0;
242   for(Int_t i=0;i<3;i++){ 
243     y[i]=fCentroid[jlin][i];
244     ey[i]=fRMSCentroid[jlin][i];
245   }
246   for(Int_t i=0;i<3;i++){
247     if(fGoodInj[jlin][i]){
248       sumY+=y[i]/ey[i]/ey[i];
249       sumX+=fPosition[i]/ey[i]/ey[i];
250       sumXX+=fPosition[i]*fPosition[i]/ey[i]/ey[i];
251       sumYY+=y[i]*y[i]/ey[i]/ey[i];
252       sumXY+=fPosition[i]*y[i]/ey[i]/ey[i];
253       sumWEI+=1./ey[i]/ey[i];
254       tzero=fTbZero/ey[i]/ey[i];
255       erry=ey[i];
256       npt++;
257     }
258   }
259   Float_t vel=0,evel=0;
260   if(npt>1){ 
261     Float_t slope=(sumWEI*sumXY-sumY*sumX)/(sumWEI*sumXX-sumX*sumX);
262     Float_t eslope=TMath::Sqrt(sumWEI/(sumWEI*sumXX-sumX*sumX));
263     vel=1./slope*10000./25.;// micron/ns
264     evel=eslope/slope/slope*10000./25.;// micron/ns
265   }
266   if(npt==1){
267     Float_t slope=(sumY-tzero)/sumX;
268     Float_t eslope=erry/sumX;
269     vel=1./slope*10000./25.;// micron/ns    
270     evel=eslope/slope/slope*10000./25.;// micron/ns
271   }
272   if(vel>fMaxDriftVel||vel<fMinDriftVel){ 
273     vel=0.;
274     evel=0.;
275   }
276   fDriftVel[jlin]=vel;
277   fSigmaDriftVel[jlin]=evel;
278 }
279 //______________________________________________________________________
280 Int_t AliITSOnlineSDDInjectors::GetAnodeNumber(Int_t iInjLine) const{
281   //
282   Int_t ian=-1;
283   if(iInjLine>32) return ian;
284   if(!fSide){
285     ian=iInjLine*8;
286     if(iInjLine==32) ian--;
287   }else{
288     ian=iInjLine*8-1;
289     if(iInjLine==0) ian=0;
290   }
291   return ian;
292 }
293
294 //______________________________________________________________________
295 void AliITSOnlineSDDInjectors::FindGoodInjectors(){
296   // 
297   for(Int_t iii=0;iii<kNInjectors;iii++){
298     Int_t ian=GetAnodeNumber(iii);
299     for(Int_t ninj=0;ninj<3;ninj++){
300       for(Int_t jjj=fTbMin[ninj];jjj<fTbMax[ninj];jjj++){
301         Float_t c1=fHisto->GetBinContent(jjj,ian+1);
302         Float_t c2=fHisto->GetBinContent(jjj+1,ian+1);
303         Float_t c3=fHisto->GetBinContent(jjj+2,ian+1);
304         if(c1>fThreshold && c2>fThreshold && c3>fThreshold){ 
305           fGoodInj[iii][ninj]=1;
306           break;
307         }
308       }
309       //      for(Int_t jjj=fTbMin[ninj];jjj<fTbMax[ninj];jjj++){
310       //        Float_t c1=fHisto->GetBinContent(jjj,ian+1);
311       //        if(c1>=fgkSaturation){
312       //          fGoodInj[iii][ninj]=0;
313       //          break;
314       //        }
315       //      }
316     }
317   }
318 }
319 //______________________________________________________________________
320 void AliITSOnlineSDDInjectors::FindCentroids(){
321   // 
322   for(Int_t iii=0;iii<kNInjectors;iii++){
323     Int_t ian=GetAnodeNumber(iii);
324     for(Int_t ninj=0;ninj<3;ninj++){
325       if(!fGoodInj[iii][ninj]) continue;
326       Float_t maxcont=0;
327       Int_t ilmax=-1;
328       for(Int_t jjj=fTbMin[ninj];jjj<fTbMax[ninj];jjj++){
329         Float_t cont=fHisto->GetBinContent(jjj,ian+1);
330         if(cont>maxcont){
331           maxcont=cont;
332           ilmax=jjj;
333         }
334       }
335       Float_t intCont=0;
336       Int_t jjj=ilmax;
337       while(1){
338         Float_t cont=fHisto->GetBinContent(jjj,ian+1);
339         if(cont<fThreshold) break;
340         if(cont<fgkSaturation){
341           fCentroid[iii][ninj]+=cont*(Float_t)jjj;
342           fRMSCentroid[iii][ninj]+=cont*TMath::Power((Float_t)jjj,2);
343           intCont+=cont;
344         }
345         jjj--;
346       }
347       jjj=ilmax+1;
348       while(1){
349         Float_t cont=fHisto->GetBinContent(jjj,ian+1);
350         if(cont<fThreshold) break;
351         if(cont<fgkSaturation){
352           fCentroid[iii][ninj]+=cont*float(jjj);
353           fRMSCentroid[iii][ninj]+=cont*TMath::Power((Float_t)jjj,2);
354           intCont+=cont;
355         }
356         jjj++;
357       }
358       if(intCont>0){ 
359         fCentroid[iii][ninj]/=intCont;
360         fRMSCentroid[iii][ninj]=TMath::Sqrt(fRMSCentroid[iii][ninj]/intCont-fCentroid[iii][ninj]*fCentroid[iii][ninj]);
361       }
362       else{ 
363         fCentroid[iii][ninj]=0.;
364         fRMSCentroid[iii][ninj]=0.;
365         fGoodInj[iii][ninj]=0;
366       }
367     }
368   }
369 }
370 //______________________________________________________________________
371 void AliITSOnlineSDDInjectors::PrintInjMap(){
372   //
373   for(Int_t iii=0;iii<kNInjectors;iii++){
374     printf("Line%d-Anode%d: %d %d %d\n",iii,GetAnodeNumber(iii),fGoodInj[iii][0],fGoodInj[iii][1],fGoodInj[iii][2]);
375   }
376 }
377 //______________________________________________________________________
378 void AliITSOnlineSDDInjectors::PrintCentroids(){
379   //
380   for(Int_t iii=0;iii<kNInjectors;iii++){
381     printf("Line%d-Anode%d: %f+-%f %f+-%f %f+-%f\n",iii,GetAnodeNumber(iii),fCentroid[iii][0],fRMSCentroid[iii][0],fCentroid[iii][1],fRMSCentroid[iii][1],fCentroid[iii][2],fRMSCentroid[iii][2]);
382   }
383 }
384 //______________________________________________________________________
385 void AliITSOnlineSDDInjectors::WriteToFXS(){
386   //
387   Char_t outfilnam[100];
388   sprintf(outfilnam,"SDDinj_mod%03d_sid%d.data",fModuleId,fSide);
389   FILE* outf=fopen(outfilnam,"w");
390   for(Int_t ic=0;ic<fPolOrder+1;ic++){
391     fprintf(outf,"%G ",fParam[ic]);
392   }
393   fprintf(outf,"\n");
394   fclose(outf);  
395 }