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