Bug fix (Andrea)
[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
4c82df4c 20/* $Id$ */
21
348f80b7 22///////////////////////////////////////////////////////////////////
23// //
24// Implementation of the class used for SDD injector analysis //
25// Origin: F.Prino, Torino, prino@to.infn.it //
26// //
27///////////////////////////////////////////////////////////////////
28
29ClassImp(AliITSOnlineSDDInjectors)
30
31const Float_t AliITSOnlineSDDInjectors::fgkSaturation = 1008.;
348f80b7 32
33//______________________________________________________________________
4ff6aa93 34AliITSOnlineSDDInjectors::AliITSOnlineSDDInjectors():AliITSOnlineSDD(),fHisto(),fTbZero(0.),fParam(),fPolOrder(0),fMinDriftVel(0.),fMaxDriftVel(0.),fThreshold(0.),fTimeDiffTB()
348f80b7 35{
36 // default constructor
37 SetMinDriftVel();
38 SetMaxDriftVel();
39 SetRangeLine1();
40 SetRangeLine2();
41 SetRangeLine3();
42 SetPositions();
43 SetPolOrder();
44 SetThreshold();
4ff6aa93 45 SetTimeDiffTB();
348f80b7 46}
47//______________________________________________________________________
979b5a5f 48AliITSOnlineSDDInjectors::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()
beb262b4 49{
50// standard constructor
348f80b7 51 SetMinDriftVel();
52 SetMaxDriftVel();
53 SetRangeLine1();
54 SetRangeLine2();
55 SetRangeLine3();
56 SetPositions();
57 SetPolOrder();
58 SetThreshold();
4ff6aa93 59 SetTimeDiffTB();
348f80b7 60}
61//______________________________________________________________________
62AliITSOnlineSDDInjectors::~AliITSOnlineSDDInjectors(){
63 // Destructor
64 if(fHisto) delete fHisto;
65 if(fParam) delete [] fParam;
66}
67//______________________________________________________________________
68void 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//______________________________________________________________________
78void AliITSOnlineSDDInjectors::Reset(){
beb262b4 79 //
348f80b7 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//______________________________________________________________________
93void 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//______________________________________________________________________
4ff6aa93 104TGraphErrors* AliITSOnlineSDDInjectors::GetLineGraph(Int_t jlin) const{
348f80b7 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//______________________________________________________________________
121Float_t AliITSOnlineSDDInjectors::GetDriftCoordinate(Float_t cAnode, Float_t cTimeBin){
beb262b4 122 //
348f80b7 123 Float_t vel=0;
124 for(Int_t i=0;i<=fPolOrder;i++) vel+=fParam[i]*TMath::Power(cAnode,(Float_t)i);
4ff6aa93 125 return vel*(cTimeBin-(fTbZero-fTimeDiffTB))*25/1000.;
348f80b7 126}
127//______________________________________________________________________
beb262b4 128TGraphErrors* AliITSOnlineSDDInjectors::GetDriftVelocityGraph() const{
348f80b7 129 //
348f80b7 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//______________________________________________________________________
142void AliITSOnlineSDDInjectors::CalcTimeBinZero(){
beb262b4 143 //
348f80b7 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//______________________________________________________________________
157void AliITSOnlineSDDInjectors::FitDriftVelocityVsAnode(){
4ff6aa93 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++){
348f80b7 164 vect[k1]=0;
4ff6aa93 165 for(Int_t k2=0;k2<kNn;k2++){
348f80b7 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 }
4ff6aa93 173 for(Int_t k1=0;k1<kNn;k1++){
348f80b7 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 }
4ff6aa93 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;
348f80b7 183 Int_t iCol=-1,iRow=-1;
4ff6aa93 184 for(Int_t i=0;i<kNn;i++){
348f80b7 185 Float_t big=0.;
4ff6aa93 186 for(Int_t j=0;j<kNn;j++){
348f80b7 187 if(iPivot[j]!=1){
4ff6aa93 188 for(Int_t k=0;k<kNn;k++){
348f80b7 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){
4ff6aa93 202 for(Int_t l=0;l<kNn;l++){
348f80b7 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;
4ff6aa93 216 for(Int_t l=0;l<kNn;l++) mat[iCol][l]*=pivinv;
348f80b7 217 vect[iCol]*=pivinv;
4ff6aa93 218 for(Int_t m=0;m<kNn;m++){
348f80b7 219 if(m!=iCol){
220 aux=mat[m][iCol];
221 mat[m][iCol]=0;
4ff6aa93 222 for(Int_t n=0;n<kNn;n++) mat[m][n]-=mat[iCol][n]*aux;
348f80b7 223 vect[m]-=vect[iCol]*aux;
224 }
225 }
226 }
227 delete [] iPivot;
228 delete [] indxR;
229 delete [] indxC;
230
231 if(fParam) delete [] fParam;
4ff6aa93 232 fParam=new Float_t[kNn];
233 for(Int_t i=0; i<kNn;i++)fParam[i]=vect[i];
348f80b7 234
4ff6aa93 235 for(Int_t i=0; i < kNn; i++) delete [] mat[i];
348f80b7 236 delete [] mat;
237 delete [] vect;
238}
239//______________________________________________________________________
240void 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//______________________________________________________________________
beb262b4 284Int_t AliITSOnlineSDDInjectors::GetAnodeNumber(Int_t iInjLine) const{
285 //
348f80b7 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}
84c6e6c0 297//______________________________________________________________________
298Int_t AliITSOnlineSDDInjectors::GetLineNumberFromAnode(Int_t nAnode) const{
299 //
300 Int_t iLine=-1;
301 if(!fSide){
302 if(nAnode%8==0) iLine=nAnode/8;
303 if(nAnode==255) iLine=32;
304 }else{
305 if(nAnode%8==7) iLine=1+nAnode/8;
306 if(nAnode==0) iLine=0;
307 }
308 if(nAnode>=256) iLine=-1;
309 return iLine;
310}
311//______________________________________________________________________
312Int_t AliITSOnlineSDDInjectors::GetAnodeStatus(Int_t nAnode) const{
313 //
314 Int_t iii=GetLineNumberFromAnode(nAnode);
315 Int_t istatus=0;
316 if(iii>=0){
317 for(Int_t ninj=0;ninj<3;ninj++) istatus+=fGoodInj[iii][ninj]<<ninj;
318 }
319 return istatus;
320}
348f80b7 321//______________________________________________________________________
322void AliITSOnlineSDDInjectors::FindGoodInjectors(){
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 for(Int_t jjj=fTbMin[ninj];jjj<fTbMax[ninj];jjj++){
328 Float_t c1=fHisto->GetBinContent(jjj,ian+1);
329 Float_t c2=fHisto->GetBinContent(jjj+1,ian+1);
330 Float_t c3=fHisto->GetBinContent(jjj+2,ian+1);
331 if(c1>fThreshold && c2>fThreshold && c3>fThreshold){
332 fGoodInj[iii][ninj]=1;
333 break;
334 }
335 }
336 // for(Int_t jjj=fTbMin[ninj];jjj<fTbMax[ninj];jjj++){
337 // Float_t c1=fHisto->GetBinContent(jjj,ian+1);
338 // if(c1>=fgkSaturation){
339 // fGoodInj[iii][ninj]=0;
340 // break;
341 // }
342 // }
343 }
344 }
345}
346//______________________________________________________________________
347void AliITSOnlineSDDInjectors::FindCentroids(){
348 //
349 for(Int_t iii=0;iii<kNInjectors;iii++){
350 Int_t ian=GetAnodeNumber(iii);
351 for(Int_t ninj=0;ninj<3;ninj++){
352 if(!fGoodInj[iii][ninj]) continue;
353 Float_t maxcont=0;
354 Int_t ilmax=-1;
355 for(Int_t jjj=fTbMin[ninj];jjj<fTbMax[ninj];jjj++){
356 Float_t cont=fHisto->GetBinContent(jjj,ian+1);
357 if(cont>maxcont){
358 maxcont=cont;
359 ilmax=jjj;
360 }
361 }
362 Float_t intCont=0;
363 Int_t jjj=ilmax;
364 while(1){
365 Float_t cont=fHisto->GetBinContent(jjj,ian+1);
366 if(cont<fThreshold) break;
367 if(cont<fgkSaturation){
368 fCentroid[iii][ninj]+=cont*(Float_t)jjj;
369 fRMSCentroid[iii][ninj]+=cont*TMath::Power((Float_t)jjj,2);
370 intCont+=cont;
371 }
372 jjj--;
373 }
374 jjj=ilmax+1;
375 while(1){
376 Float_t cont=fHisto->GetBinContent(jjj,ian+1);
377 if(cont<fThreshold) break;
378 if(cont<fgkSaturation){
379 fCentroid[iii][ninj]+=cont*float(jjj);
380 fRMSCentroid[iii][ninj]+=cont*TMath::Power((Float_t)jjj,2);
381 intCont+=cont;
382 }
383 jjj++;
384 }
385 if(intCont>0){
386 fCentroid[iii][ninj]/=intCont;
387 fRMSCentroid[iii][ninj]=TMath::Sqrt(fRMSCentroid[iii][ninj]/intCont-fCentroid[iii][ninj]*fCentroid[iii][ninj]);
388 }
389 else{
390 fCentroid[iii][ninj]=0.;
391 fRMSCentroid[iii][ninj]=0.;
392 fGoodInj[iii][ninj]=0;
393 }
4c82df4c 394 if(fRMSCentroid[iii][ninj]==0) fGoodInj[iii][ninj]=0;
348f80b7 395 }
396 }
397}
398//______________________________________________________________________
399void AliITSOnlineSDDInjectors::PrintInjMap(){
400 //
401 for(Int_t iii=0;iii<kNInjectors;iii++){
402 printf("Line%d-Anode%d: %d %d %d\n",iii,GetAnodeNumber(iii),fGoodInj[iii][0],fGoodInj[iii][1],fGoodInj[iii][2]);
403 }
404}
405//______________________________________________________________________
406void AliITSOnlineSDDInjectors::PrintCentroids(){
407 //
408 for(Int_t iii=0;iii<kNInjectors;iii++){
409 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]);
410 }
411}
412//______________________________________________________________________
4c82df4c 413void AliITSOnlineSDDInjectors::WriteToASCII(Int_t evNumb, UInt_t timeStamp, Int_t optAppend){
348f80b7 414 //
415 Char_t outfilnam[100];
979b5a5f 416 sprintf(outfilnam,"SDDinj_ddl%02dc%02d_sid%d.data",fDDL,fCarlos,fSide);
4c82df4c 417 FILE* outf;
18da6e54 418 if(optAppend==0){
419 outf=fopen(outfilnam,"w");
420 fprintf(outf,"%d\n",fPolOrder);
421 }
4c82df4c 422 else outf=fopen(outfilnam,"a");
423 fprintf(outf,"%d %d ",evNumb,timeStamp);
348f80b7 424 for(Int_t ic=0;ic<fPolOrder+1;ic++){
425 fprintf(outf,"%G ",fParam[ic]);
426 }
427 fprintf(outf,"\n");
428 fclose(outf);
429}