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