First version of the SDD DA calibration classes. AliITSOnlineSDDBase - for measuremen...
[u/mrichter/AliRoot.git] / ITS / AliITSOnlineSDDInjectors.cxx
CommitLineData
348f80b7 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
27ClassImp(AliITSOnlineSDDInjectors)
28
29const Float_t AliITSOnlineSDDInjectors::fgkSaturation = 1008.;
30const 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//______________________________________________________________________
46AliITSOnlineSDDInjectors::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//______________________________________________________________________
58AliITSOnlineSDDInjectors::~AliITSOnlineSDDInjectors(){
59 // Destructor
60 if(fHisto) delete fHisto;
61 if(fParam) delete [] fParam;
62}
63//______________________________________________________________________
64void 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//______________________________________________________________________
74void 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//______________________________________________________________________
88void 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//______________________________________________________________________
99TGraphErrors* 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//______________________________________________________________________
116Float_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//______________________________________________________________________
122TGraphErrors* 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//______________________________________________________________________
137void 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//______________________________________________________________________
151void 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//______________________________________________________________________
233void 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//______________________________________________________________________
277Int_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//______________________________________________________________________
291void 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//______________________________________________________________________
316void 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//______________________________________________________________________
367void 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//______________________________________________________________________
374void 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//______________________________________________________________________
381void 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}