SPD services on the cones (A. Pulvirenti)
[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
19f1458a 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;
348f80b7 67 if(fParam) delete [] fParam;
68}
69//______________________________________________________________________
70void 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//______________________________________________________________________
80void AliITSOnlineSDDInjectors::Reset(){
beb262b4 81 //
348f80b7 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//______________________________________________________________________
95void 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//______________________________________________________________________
4ff6aa93 106TGraphErrors* AliITSOnlineSDDInjectors::GetLineGraph(Int_t jlin) const{
348f80b7 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//______________________________________________________________________
123Float_t AliITSOnlineSDDInjectors::GetDriftCoordinate(Float_t cAnode, Float_t cTimeBin){
beb262b4 124 //
348f80b7 125 Float_t vel=0;
126 for(Int_t i=0;i<=fPolOrder;i++) vel+=fParam[i]*TMath::Power(cAnode,(Float_t)i);
4ff6aa93 127 return vel*(cTimeBin-(fTbZero-fTimeDiffTB))*25/1000.;
348f80b7 128}
129//______________________________________________________________________
beb262b4 130TGraphErrors* AliITSOnlineSDDInjectors::GetDriftVelocityGraph() const{
348f80b7 131 //
348f80b7 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//______________________________________________________________________
144void AliITSOnlineSDDInjectors::CalcTimeBinZero(){
beb262b4 145 //
348f80b7 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//______________________________________________________________________
159void AliITSOnlineSDDInjectors::FitDriftVelocityVsAnode(){
4ff6aa93 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++){
348f80b7 166 vect[k1]=0;
4ff6aa93 167 for(Int_t k2=0;k2<kNn;k2++){
348f80b7 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 }
4ff6aa93 175 for(Int_t k1=0;k1<kNn;k1++){
348f80b7 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 }
4ff6aa93 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;
348f80b7 185 Int_t iCol=-1,iRow=-1;
4ff6aa93 186 for(Int_t i=0;i<kNn;i++){
348f80b7 187 Float_t big=0.;
4ff6aa93 188 for(Int_t j=0;j<kNn;j++){
348f80b7 189 if(iPivot[j]!=1){
4ff6aa93 190 for(Int_t k=0;k<kNn;k++){
348f80b7 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){
4ff6aa93 204 for(Int_t l=0;l<kNn;l++){
348f80b7 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;
4ff6aa93 218 for(Int_t l=0;l<kNn;l++) mat[iCol][l]*=pivinv;
348f80b7 219 vect[iCol]*=pivinv;
4ff6aa93 220 for(Int_t m=0;m<kNn;m++){
348f80b7 221 if(m!=iCol){
222 aux=mat[m][iCol];
223 mat[m][iCol]=0;
4ff6aa93 224 for(Int_t n=0;n<kNn;n++) mat[m][n]-=mat[iCol][n]*aux;
348f80b7 225 vect[m]-=vect[iCol]*aux;
226 }
227 }
228 }
229 delete [] iPivot;
230 delete [] indxR;
231 delete [] indxC;
232
233 if(fParam) delete [] fParam;
4ff6aa93 234 fParam=new Float_t[kNn];
235 for(Int_t i=0; i<kNn;i++)fParam[i]=vect[i];
348f80b7 236
4ff6aa93 237 for(Int_t i=0; i < kNn; i++) delete [] mat[i];
348f80b7 238 delete [] mat;
239 delete [] vect;
240}
241//______________________________________________________________________
242void 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++){
53855521 253 if(fGoodInj[jlin][i] && ey[i]!=0){
348f80b7 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));
53855521 269 if(slope!=0){
270 vel=1./slope*10000./25.;// micron/ns
271 evel=eslope/slope/slope*10000./25.;// micron/ns
272 }
348f80b7 273 }
274 if(npt==1){
275 Float_t slope=(sumY-tzero)/sumX;
276 Float_t eslope=erry/sumX;
53855521 277 if(slope!=0){
278 vel=1./slope*10000./25.;// micron/ns
279 evel=eslope/slope/slope*10000./25.;// micron/ns
280 }
348f80b7 281 }
282 if(vel>fMaxDriftVel||vel<fMinDriftVel){
283 vel=0.;
284 evel=0.;
285 }
286 fDriftVel[jlin]=vel;
287 fSigmaDriftVel[jlin]=evel;
288}
289//______________________________________________________________________
beb262b4 290Int_t AliITSOnlineSDDInjectors::GetAnodeNumber(Int_t iInjLine) const{
291 //
348f80b7 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}
84c6e6c0 303//______________________________________________________________________
304Int_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//______________________________________________________________________
318Int_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}
348f80b7 327//______________________________________________________________________
328void 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//______________________________________________________________________
353void 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 }
4c82df4c 400 if(fRMSCentroid[iii][ninj]==0) fGoodInj[iii][ninj]=0;
348f80b7 401 }
402 }
403}
404//______________________________________________________________________
405void 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//______________________________________________________________________
412void 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//______________________________________________________________________
4c82df4c 419void AliITSOnlineSDDInjectors::WriteToASCII(Int_t evNumb, UInt_t timeStamp, Int_t optAppend){
348f80b7 420 //
421 Char_t outfilnam[100];
979b5a5f 422 sprintf(outfilnam,"SDDinj_ddl%02dc%02d_sid%d.data",fDDL,fCarlos,fSide);
4c82df4c 423 FILE* outf;
18da6e54 424 if(optAppend==0){
425 outf=fopen(outfilnam,"w");
426 fprintf(outf,"%d\n",fPolOrder);
427 }
4c82df4c 428 else outf=fopen(outfilnam,"a");
429 fprintf(outf,"%d %d ",evNumb,timeStamp);
348f80b7 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}