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