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