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 <TFile.h>
  16 #include "AliITSOnlineSDDInjectors.h"
  17 #include "AliLog.h"
  18 #include <TH2F.h>
  19 #include <TF1.h>
  20 #include <TGraphErrors.h>
  21 #include <TMath.h>
  22 #include <TString.h>
  23
  24 /* $Id$ */
  25
  26 ///////////////////////////////////////////////////////////////////
  27 //                                                               //
  28 // Implementation of the class used for SDD injector analysis    //
  29 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  30 //                                                               //
  31 ///////////////////////////////////////////////////////////////////
  32
  33 ClassImp(AliITSOnlineSDDInjectors)
  34
  35 const Float_t AliITSOnlineSDDInjectors::fgkSaturation = 1008.;
  36 const Float_t AliITSOnlineSDDInjectors::fgkDefaultLThreshold1 = 8.;
  37 const Float_t AliITSOnlineSDDInjectors::fgkDefaultLThreshold = 15.;
  38 const Float_t AliITSOnlineSDDInjectors::fgkDefaultHThreshold1 =15.;
  39 const Float_t AliITSOnlineSDDInjectors::fgkDefaultHThreshold = 30.;
  40 const Float_t AliITSOnlineSDDInjectors::fgkDefaultMinSpeed = 5.5;
  41 const Float_t AliITSOnlineSDDInjectors::fgkDefaultMaxSpeed = 9.0;
  42 const Float_t AliITSOnlineSDDInjectors::fgkDefaultMaxErr = 1.5;
  43 const Int_t   AliITSOnlineSDDInjectors::fgkDefaultPolDegree = 3;
  44 const Float_t AliITSOnlineSDDInjectors::fgkDefaultTimeStep = 50.;
  45 const UShort_t AliITSOnlineSDDInjectors::fgkDefaultTbMin[kInjLines] = {10,50,100};
  46 const UShort_t AliITSOnlineSDDInjectors::fgkDefaultTbMax[kInjLines] = {20,70,120};
  47
  48 //______________________________________________________________________
  49 AliITSOnlineSDDInjectors::AliITSOnlineSDDInjectors():
  50   AliITSOnlineSDD(),
  51   fHisto(),
  52   fTbZero(0.),
  53   fRMSTbZero(0.),
  54   fNEvents(0),
  55   fParam(),
  56   fPolDegree(0),
  57   fActualPolDegree(0),
  58   fMinDriftSpeed(0.),
  59   fMaxDriftSpeed(0.),
  60   fMaxDriftSpeedErr(0.),
  61   fFirstPadForFit(0),
  62   fLastPadForFit(0),
  63   fPadStatusCutForFit(0),
  64   fTimeStep(0.),
  65   fUseTimeZeroSignal(kFALSE),
  66   fMaxCellsAboveThreshold(40)
  67 {
  68   // default constructor
  69   SetPositions();
  70   SetDefaults();
  71   SetTimeStep(fgkDefaultTimeStep);
  72   for(Int_t i=0;i<kInjPads;i++){
  73     fSumDriftSpeed[i]=0.;
  74     fSumSqDriftSpeed[i]=0.;
  75     fSumPadStatus[i]=0;
  76     fSumPadStatusCut[i]=0;
  77     fNEventsInPad[i]=0;
  78   }
  79   Reset();
  80 }
  81 //______________________________________________________________________
  82 AliITSOnlineSDDInjectors::AliITSOnlineSDDInjectors(Int_t nddl, Int_t ncarlos, Int_t sid):
  83   AliITSOnlineSDD(nddl,ncarlos,sid),
  84   fHisto(),
  85   fTbZero(0.),
  86   fRMSTbZero(0.),
  87   fNEvents(0),
  88   fParam(),
  89   fPolDegree(0),
  90   fActualPolDegree(0),
  91   fMinDriftSpeed(0.),
  92   fMaxDriftSpeed(0.),
  93   fMaxDriftSpeedErr(0.),
  94   fFirstPadForFit(0),
  95   fLastPadForFit(0),
  96   fPadStatusCutForFit(0),
  97   fTimeStep(0.),
  98   fUseTimeZeroSignal(kFALSE),
  99   fMaxCellsAboveThreshold(40)
 100 {
 101 // standard constructor
 102   SetPositions();
 103   SetDefaults();
 104   SetTimeStep(fgkDefaultTimeStep);
 105   for(Int_t i=0;i<kInjPads;i++){
 106     fSumDriftSpeed[i]=0.;
 107     fSumSqDriftSpeed[i]=0.;
 108     fSumPadStatus[i]=0;
 109     fSumPadStatusCut[i]=0;
 110     fNEventsInPad[i]=0;
 111   }
 112   Reset();
 113 }
 114 //______________________________________________________________________
 115 AliITSOnlineSDDInjectors::~AliITSOnlineSDDInjectors(){
 116   // Destructor
 117   // fHisto should not be deleted here because it points to an histo created
 118   // by the external code which calls the method AnalyzeEvent
 119   // if(fHisto) delete fHisto;
 120   if(fParam) delete [] fParam;
 121 }
 122 //______________________________________________________________________
 123 void AliITSOnlineSDDInjectors::SetDefaults(){
 124   // Sets default values for parameters
 125   for(Int_t i=0;i<kInjLines;i++) {
 126     SetInjLineRange(i,fgkDefaultTbMin[i],fgkDefaultTbMax[i]);
 127     SetUseLine(i,kTRUE);
 128     SetThresholds(i,fgkDefaultLThreshold,fgkDefaultHThreshold);
 129   }
 130   SetThresholds(0,fgkDefaultLThreshold1,fgkDefaultHThreshold1);
 131   SetPolDegree(fgkDefaultPolDegree);
 132   SetMinDriftSpeed(fgkDefaultMinSpeed);
 133   SetMaxDriftSpeed(fgkDefaultMaxSpeed);
 134   SetMaxDriftSpeedErr(fgkDefaultMaxErr);
 135   SetFitLimits(1,kInjPads-2); // exclude first and last pad
 136   SetPadStatusCutForFit();
 137 }
 138 //______________________________________________________________________
 139 void AliITSOnlineSDDInjectors::Set20MHzConfig(){
 140   // Sets specific parameters for 20 MHz running
 141   SetInjLineRange(0,10,20);
 142   SetInjLineRange(1,50,70);
 143   SetInjLineRange(2,100,120);
 144   SetTimeStep(50.);
 145   SetMaxNumberOfCellsPerAnode(40);
 146 }
 147 //______________________________________________________________________
 148 void AliITSOnlineSDDInjectors::Set40MHzConfig(){
 149   // Sets specific parameters for 20 MHz running
 150   SetInjLineRange(0,20,50);
 151   SetInjLineRange(1,90,160);
 152   SetInjLineRange(2,170,240);
 153   SetTimeStep(25.);
 154   SetMaxNumberOfCellsPerAnode(80);
 155 }
 156 //______________________________________________________________________
 157 void AliITSOnlineSDDInjectors::SetPositions(){
 158   // Sets drift distances for the 3 injector lines
 159   Double_t xLinFromCenterUm[kInjLines]={31860.,17460.,660.};
 160   Double_t xAnodeFromCenterUm=35085;
 161   for(Int_t i=0;i<kInjLines;i++){
 162     fPosition[i]=xAnodeFromCenterUm-xLinFromCenterUm[i];
 163     fPosition[i]/=10000.; // from microns to cm
 164   }
 165 }
 166 //______________________________________________________________________
 167 void AliITSOnlineSDDInjectors::Reset(){
 168   // Resets all counters
 169   for(Int_t i=0;i<kInjPads;i++){
 170     fDriftSpeed[i]=0.;
 171     fDriftSpeedErr[i]=0.;
 172   }
 173   for(Int_t i=0;i<kInjPads;i++){
 174     for(Int_t j=0;j<kInjLines;j++){
 175       fGoodInj[i][j]=0;
 176       fCentroid[i][j]=0.;
 177       fRMSCentroid[i][j]=0.;
 178     }
 179   }
 180 }
 181 //______________________________________________________________________
 182 void AliITSOnlineSDDInjectors::AnalyzeEvent(TH2F* his){
 183   // Analyze the current event
 184   AddEvent(his);
 185   FitDriftSpeedVsAnode();
 186 }
 187 //______________________________________________________________________
 188 void AliITSOnlineSDDInjectors::AddEvent(TH2F* his){
 189   // Add the drift speed from current event to the average value
 190   if(fNEvents==0){
 191     for(Int_t i=0;i<kInjPads;i++){
 192       fSumDriftSpeed[i]=0.;
 193       fSumSqDriftSpeed[i]=0.;
 194       fSumPadStatus[i]=0;
 195       fSumPadStatusCut[i]=0;
 196       fNEventsInPad[i]=0;
 197     }
 198   }
 199   Reset();
 200   fHisto=his;
 201   FindGoodInjectors();
 202   FindCentroids();
 203   CalcTimeBinZero();
 204   for(Int_t j=0;j<kInjPads;j++){
 205     CalcDriftSpeed(j);
 206     Int_t padStatus=GetInjPadStatus(j);
 207     fSumPadStatus[j]+=padStatus;
 208     if(padStatus>fPadStatusCutForFit){
 209       fSumDriftSpeed[j]+=fDriftSpeed[j];
 210       fSumSqDriftSpeed[j]+=fDriftSpeed[j]*fDriftSpeed[j];
 211       fSumPadStatusCut[j]+=padStatus;
 212       fNEventsInPad[j]++;
 213     }
 214   }
 215   ++fNEvents;
 216 }
 217 //______________________________________________________________________
 218 Double_t AliITSOnlineSDDInjectors::GetRMSDriftSpeed(Int_t ipad) const {
 219   // Compute RMS of drift speed distribution on one anode
 220   if(fNEventsInPad[ipad]<=1) return 0.;
 221   Double_t mean=fSumDriftSpeed[ipad]/(Double_t)fNEventsInPad[ipad];
 222   Double_t diff=fSumSqDriftSpeed[ipad]/(Double_t)fNEventsInPad[ipad]-mean*mean;
 223   if(diff<0.) diff=0.;
 224   return TMath::Sqrt(diff);
 225 }
 226
 227 //______________________________________________________________________
 228 void AliITSOnlineSDDInjectors::FitMeanDriftSpeedVsAnode(){
 229   // Fits the average drift speed vs.anode number
 230   if(fNEvents==0) return;
 231   for(Int_t i=0;i<kInjPads;i++){
 232     fDriftSpeed[i]=GetMeanDriftSpeed(i);
 233     Int_t padStatusCut=(Int_t)(GetMeanPadStatusCut(i)+0.5);
 234     for(Int_t ilin=0; ilin<kInjLines ; ilin++) fGoodInj[i][ilin]=(padStatusCut&1<<ilin)>>ilin;
 235     if(fNEventsInPad[i]>1){
 236       Double_t rms=GetRMSDriftSpeed(i);
 237       if(rms>0.) fDriftSpeedErr[i]=rms/TMath::Sqrt(fNEventsInPad[i]);
 238     }else{
 239       for(Int_t ilin=0; ilin<kInjLines ; ilin++) fGoodInj[i][ilin]=0;
 240     }
 241   }
 242   FitDriftSpeedVsAnode();
 243   for(Int_t i=0;i<kInjPads;i++){
 244     Int_t padStatus=(Int_t)(GetMeanPadStatusCut(i)+0.5);
 245     for(Int_t ilin=0; ilin<kInjLines ; ilin++) fGoodInj[i][ilin]=(padStatus&1<<ilin)>>ilin;
 246   }
 247 }
 248 //______________________________________________________________________
 249 TGraphErrors* AliITSOnlineSDDInjectors::GetTimeVsDistGraph(Int_t jpad) const{
 250   // Builds the graph of drift time vs. drift distance
 251   const Int_t kPts=kInjLines+1;
 252   Float_t x[kPts],y[kPts],ex[kPts],ey[kPts];
 253   x[0]=0.;
 254   ex[0]=0.;
 255   y[0]=fTbZero;
 256   ey[0]=0.;
 257   for(Int_t i=0;i<kInjLines;i++){
 258     x[i+1]=fPosition[i];
 259     ex[i+1]=0.;
 260     y[i+1]=fCentroid[jpad][i];
 261     ey[i+1]=fRMSCentroid[jpad][i];
 262   }
 263   TGraphErrors *g=new TGraphErrors(4,x,y,ex,ey);
 264   return g;
 265 }
 266
 267 //______________________________________________________________________
 268 TGraphErrors* AliITSOnlineSDDInjectors::GetDriftSpeedGraph() const{
 269   // Builds the graph of drift speed vs. anode number
 270   Int_t ipt=0;
 271   TGraphErrors *g=new TGraphErrors(0);
 272   for(Int_t i=0;i<kInjPads;i++){
 273     if(fDriftSpeed[i]>0){
 274       g->SetPoint(ipt,GetAnodeNumber(i),fDriftSpeed[i]);
 275       g->SetPointError(ipt,0,fDriftSpeedErr[i]);
 276       ipt++;
 277     }
 278   }
 279   return g;
 280 }
 281 //______________________________________________________________________
 282 TGraphErrors* AliITSOnlineSDDInjectors::GetSelectedDriftSpeedGraph(Int_t minAcceptStatus) const{
 283   // TGraphErrors with only pads with status of injector >= minAcceptStatus
 284   Int_t ipt=0;
 285   TGraphErrors *g=new TGraphErrors(0);
 286   for(Int_t i=0;i<kInjPads;i++){
 287     Int_t padStatus = GetInjPadStatus(i);
 288     if(fDriftSpeed[i]>0 && padStatus >= minAcceptStatus ){
 289       g->SetPoint(ipt,GetAnodeNumber(i),fDriftSpeed[i]);
 290       g->SetPointError(ipt,0,fDriftSpeedErr[i]);
 291       ipt++;
 292     }
 293   }
 294   return g;
 295 }
 296 //______________________________________________________________________
 297 void AliITSOnlineSDDInjectors::CalcTimeBinZero(){
 298   // Get time zero from trigger signal
 299   Double_t tzero=0.,intCont=0.,rmsPeak=0.;
 300   Bool_t isTbUsed[256];
 301   Int_t nTbUsed=0;
 302   for(Int_t i=0;i<256;i++) isTbUsed[i]=0;
 303   for(Int_t ian=0;ian<fgkNAnodes;ian++){
 304     for(Int_t itb=1;itb<fTbMin[0];itb++){
 305       Double_t cont=fHisto->GetBinContent(itb,ian+1);
 306       Double_t contm1=fHisto->GetBinContent(itb+1,ian+1);
 307       Double_t contp1=fHisto->GetBinContent(itb-1,ian+1);
 308       if(cont>fLowThreshold[0]){
 309         if(cont>fHighThreshold[0] &&(contm1>fLowThreshold[0] || contp1>fLowThreshold[0])){
 310           tzero+=cont*float(itb);
 311           rmsPeak+=cont*float(itb)*float(itb);
 312           intCont+=cont;
 313           if(!isTbUsed[itb]){
 314             isTbUsed[itb]=1;
 315             ++nTbUsed;
 316           }
 317         }
 318       }
 319     }
 320   }
 321   if(intCont>0){
 322     fTbZero=tzero/intCont;
 323     fRMSTbZero=TMath::Sqrt(rmsPeak/intCont-fTbZero*fTbZero);
 324   }
 325   if(nTbUsed==1) fRMSTbZero=0.5;
 326 }
 327
 328 //______________________________________________________________________
 329 void AliITSOnlineSDDInjectors::FitDriftSpeedVsAnode(){
 330   // fits the anode dependence of drift speed
 331
 332   Float_t rangeForMax[2]={78.,178.};
 333   PolyFit(fPolDegree);
 334   fActualPolDegree=fPolDegree;
 335   if(fPolDegree==3){
 336     Double_t deltasq=fParam[2]*fParam[2]-3*fParam[1]*fParam[3];
 337     Double_t zero1=-999.;
 338     Double_t zero2=-999.;
 339     if(deltasq>=0. && TMath::Abs(fParam[3])>0.){
 340       Double_t delta=TMath::Sqrt(deltasq);
 341       zero1=(-fParam[2]+delta)/3./fParam[3];
 342       zero2=(-fParam[2]-delta)/3./fParam[3];
 343     }
 344     Bool_t twoZeroes=kFALSE;
 345     Bool_t oneZero=kFALSE;
 346     if(zero1>0. && zero1<256. && zero2>0. && zero2<256.) twoZeroes=kTRUE;
 347     if(zero1>rangeForMax[0] && zero1<rangeForMax[1]) oneZero=kTRUE;
 348     if(zero2>rangeForMax[0] && zero2<rangeForMax[1]) oneZero=kTRUE;
 349     if(!oneZero || twoZeroes){
 350       PolyFit(2);
 351       Double_t xmax=-999.;
 352       if(fParam[2]<0.) xmax=-fParam[1]/2./fParam[2];
 353       if(xmax>rangeForMax[0] && xmax<rangeForMax[1]){
 354         fActualPolDegree=2;
 355       }else{
 356         Double_t averSpeed=0.;
 357         Double_t sumWei=0.;
 358         Int_t nUsedPts=0;
 359         for(Int_t jpad=fFirstPadForFit; jpad<=fLastPadForFit; jpad++){
 360           if(fDriftSpeed[jpad]>0 && GetInjPadStatus(jpad)>fPadStatusCutForFit){
 361             Double_t wei=1./fDriftSpeedErr[jpad]/fDriftSpeedErr[jpad];
 362             averSpeed+=wei*fDriftSpeed[jpad];
 363             sumWei+=wei;
 364             nUsedPts++;
 365           }
 366         }
 367         if(sumWei>0.) averSpeed/=sumWei;
 368         if(nUsedPts<fPolDegree+1) averSpeed=0;
 369         fParam[0]=averSpeed;
 370         for(Int_t i=1; i < fPolDegree+1; i++) fParam[i]=0.;
 371         fActualPolDegree=0;
 372       }
 373     }
 374   }
 375 }
 376 //______________________________________________________________________
 377 void AliITSOnlineSDDInjectors::PolyFit(Int_t degree){
 378   // fits the anode dependence of drift speed with a polynomial function
 379   const Int_t kNn=degree+1;
 380   const Int_t kDimens=fPolDegree+1;
 381
 382   Double_t **mat = new Double_t*[kNn];
 383   for(Int_t i=0; i < kNn; i++) mat[i] = new Double_t[kNn];
 384   Double_t *vect = new Double_t[kNn];
 385
 386   for(Int_t k1=0;k1<kNn;k1++){
 387     vect[k1]=0;
 388     for(Int_t k2=0;k2<kNn;k2++){
 389       mat[k1][k2]=0;
 390     }
 391   }
 392   Int_t npts = 0;
 393   for(Int_t k1=0;k1<kNn;k1++){
 394     for(Int_t jpad=fFirstPadForFit; jpad<=fLastPadForFit; jpad++){
 395       Double_t x=(Double_t)GetAnodeNumber(jpad);
 396       if(fDriftSpeed[jpad]>0 && GetInjPadStatus(jpad)>fPadStatusCutForFit){
 397         vect[k1]+=fDriftSpeed[jpad]*TMath::Power(x,k1)/TMath::Power(fDriftSpeedErr[jpad],2);
 398         if(k1==0) npts++;
 399         for(Int_t k2=0;k2<kNn;k2++){
 400           mat[k1][k2]+=TMath::Power(x,k1+k2)/TMath::Power(fDriftSpeedErr[jpad],2);
 401         }
 402       }
 403     }
 404   }
 405   if(npts<fPolDegree+1){
 406     if(fParam) delete [] fParam;
 407     fParam=new Double_t[kDimens];
 408     for(Int_t i=0; i<kDimens;i++)fParam[i]=0;
 409   }else{
 410     Int_t *iPivot = new Int_t[kNn];
 411     Int_t *indxR = new Int_t[kNn];
 412     Int_t *indxC = new Int_t[kNn];
 413     for(Int_t i=0;i<kNn;i++) iPivot[i]=0;
 414     Int_t iCol=-1,iRow=-1;
 415     for(Int_t i=0;i<kNn;i++){
 416       Double_t big=0.;
 417       for(Int_t j=0;j<kNn;j++){
 418         if(iPivot[j]!=1){
 419           for(Int_t k=0;k<kNn;k++){
 420             if(iPivot[k]==0){
 421               if(TMath::Abs(mat[j][k])>=big){
 422                 big=TMath::Abs(mat[j][k]);
 423                 iRow=j;
 424                 iCol=k;
 425               }
 426             }
 427           }
 428         }
 429       }
 430       iPivot[iCol]++;
 431       Double_t aux;
 432       if(iRow!=iCol){
 433         for(Int_t l=0;l<kNn;l++){
 434           aux=mat[iRow][l];
 435           mat[iRow][l]=mat[iCol][l];
 436           mat[iCol][l]=aux;
 437         }
 438         aux=vect[iRow];
 439         vect[iRow]=vect[iCol];
 440         vect[iCol]=aux;
 441       }
 442       indxR[i]=iRow;
 443       indxC[i]=iCol;
 444       if(mat[iCol][iCol]==0) break;
 445       Double_t pivinv=1./mat[iCol][iCol];
 446       mat[iCol][iCol]=1;
 447       for(Int_t l=0;l<kNn;l++) mat[iCol][l]*=pivinv;
 448       vect[iCol]*=pivinv;
 449       for(Int_t m=0;m<kNn;m++){
 450         if(m!=iCol){
 451           aux=mat[m][iCol];
 452           mat[m][iCol]=0;
 453           for(Int_t n=0;n<kNn;n++) mat[m][n]-=mat[iCol][n]*aux;
 454           vect[m]-=vect[iCol]*aux;
 455         }
 456       }
 457     }
 458     delete [] iPivot;
 459     delete [] indxR;
 460     delete [] indxC;
 461
 462
 463     if(fParam) delete [] fParam;
 464     fParam=new Double_t[kDimens];
 465     for(Int_t i=0; i<kNn;i++)fParam[i]=vect[i];
 466     if(degree<fPolDegree) for(Int_t i=kNn; i<kDimens;i++)fParam[i]=0.;
 467   }
 468
 469   for(Int_t i=0; i < kNn; i++) delete [] mat[i];
 470   delete [] mat;
 471   delete [] vect;
 472 }
 473 //______________________________________________________________________
 474 void AliITSOnlineSDDInjectors::CalcDriftSpeed(Int_t jpad){
 475   // Computes the drift speed from the fit to the 3 injector lines for each anode
 476   Double_t sumY=0,sumX=0,sumXX=0,sumYY=0.,sumXY=0,sumWEI=0.;
 477   Int_t npt=0;
 478   Double_t y[kInjLines],ey[kInjLines];
 479   Double_t tzero=0,erry=0;
 480   for(Int_t i=0;i<kInjLines;i++){
 481     y[i]=fCentroid[jpad][i];
 482     ey[i]=fRMSCentroid[jpad][i];
 483   }
 484   for(Int_t i=0;i<kInjLines;i++){
 485     if(!fUseLine[i]) continue;
 486     if(fGoodInj[jpad][i] && ey[i]!=0){
 487       sumY+=y[i]/ey[i]/ey[i];
 488       sumX+=fPosition[i]/ey[i]/ey[i];
 489       sumXX+=fPosition[i]*fPosition[i]/ey[i]/ey[i];
 490       sumYY+=y[i]*y[i]/ey[i]/ey[i];
 491       sumXY+=fPosition[i]*y[i]/ey[i]/ey[i];
 492       sumWEI+=1./ey[i]/ey[i];
 493       tzero=fTbZero/ey[i]/ey[i];
 494       erry=ey[i]/ey[i]/ey[i];
 495       npt++;
 496     }
 497   }
 498   Double_t slope=0.,eslope=0.;
 499   if(npt==1){
 500     slope=(sumY-tzero)/sumX;
 501     eslope=erry/sumX;
 502   }
 503   if(npt>1){
 504     if(fUseTimeZeroSignal){
 505       sumY+=fTbZero/fRMSTbZero/fRMSTbZero;
 506       sumX+=0.;
 507       sumXX+=0.;
 508       sumYY+=fTbZero*fTbZero/fRMSTbZero/fRMSTbZero;
 509       sumXY+=0.;
 510       sumWEI+=1./fRMSTbZero/fRMSTbZero;
 511     }
 512     slope=(sumWEI*sumXY-sumY*sumX)/(sumWEI*sumXX-sumX*sumX);
 513     eslope=TMath::Sqrt(sumWEI/(sumWEI*sumXX-sumX*sumX));
 514   }
 515
 516   Double_t vel=0,evel=0;
 517   if(slope!=0. && fTimeStep>0.){
 518     vel=1./slope*10000./fTimeStep;// micron/ns
 519     evel=eslope/slope/slope*10000./fTimeStep;// micron/ns
 520   }
 521   if(vel>fMaxDriftSpeed||vel<fMinDriftSpeed || evel>fMaxDriftSpeedErr){
 522     vel=0.;
 523     evel=0.;
 524   }
 525   fDriftSpeed[jpad]=vel;
 526   fDriftSpeedErr[jpad]=evel;
 527 }
 528 //______________________________________________________________________
 529 Int_t AliITSOnlineSDDInjectors::GetAnodeNumber(Int_t iInjPad) const{
 530   // Injectors location along anodes:
 531   // Side left  (UP)   - channel 0: injectors on anodes 0,7,15,...,247,255
 532   // Side right (DOWN) - channel 1: injectors on anodes 0,8,16,...,248,255
 533   Int_t ian=-1;
 534   if(iInjPad>=kInjPads) return ian;
 535   if(fSide==1){  // right side
 536     ian=iInjPad*8;
 537     if(iInjPad==32) ian--;
 538   }else{         // left side
 539     ian=iInjPad*8-1;
 540     if(iInjPad==0) ian=0;
 541   }
 542   return ian;
 543 }
 544 //______________________________________________________________________
 545 Int_t AliITSOnlineSDDInjectors::GetInjPadNumberFromAnode(Int_t nAnode) const{
 546   // Converts anode number into injector pad index
 547   Int_t iInjPad=-1;
 548   if(fSide==1){  // right side
 549     if(nAnode%8==0) iInjPad=nAnode/8;
 550     if(nAnode==255) iInjPad=32;
 551   }else{         // left side
 552     if(nAnode%8==7) iInjPad=1+nAnode/8;
 553     if(nAnode==0) iInjPad=0;
 554   }
 555   if(nAnode>=256) iInjPad=-1;
 556   return iInjPad;
 557 }
 558 //______________________________________________________________________
 559 Int_t AliITSOnlineSDDInjectors::GetInjPadStatus(Int_t jpad) const{
 560   // returns an integer value with status of injector lines for given pad/anode
 561   // status=7  -->  111  all injector are good
 562   // status=6  -->  110  1st line (close to anodes) is bad, other two are good
 563   // ....
 564   // status=1  -->  001  only 1st line (close to anodes) good
 565   // status=0  -->  000  all lines are bad
 566   Int_t istatus=0;
 567   if(jpad>=0 && jpad<kInjPads){
 568     for(Int_t jlin=0;jlin<kInjLines;jlin++) istatus+=fGoodInj[jpad][jlin]<<jlin;
 569   }
 570   return istatus;
 571 }
 572 //______________________________________________________________________
 573 void AliITSOnlineSDDInjectors::FindGoodInjectors(){
 574   // Mark good injector pads
 575   // good = 1 cell above high threshold + 1 neighbour above low threshold
 576   for(Int_t jpad=0;jpad<kInjPads;jpad++){
 577     Int_t ian=GetAnodeNumber(jpad);
 578     Int_t countAbove=0;
 579     for(Int_t jjj=0; jjj<fHisto->GetNbinsX(); jjj++){
 580       Float_t c=fHisto->GetBinContent(jjj+1,ian+1);
 581       if(c>0.5) countAbove++;
 582     }
 583     if(countAbove>fMaxCellsAboveThreshold){
 584       for(Int_t jlin=0;jlin<kInjLines;jlin++) fGoodInj[jpad][jlin]=0;
 585     }else{
 586       for(Int_t jlin=0;jlin<kInjLines;jlin++){
 587         for(Int_t jjj=fTbMin[jlin];jjj<fTbMax[jlin];jjj++){
 588           Float_t c1=fHisto->GetBinContent(jjj,ian+1);
 589           Float_t c2=fHisto->GetBinContent(jjj+1,ian+1);
 590           //    Float_t c3=fHisto->GetBinContent(jjj+2,ian+1);
 591           if(c1>fLowThreshold[jlin] && c2>fLowThreshold[jlin]){
 592             if(c1>fHighThreshold[jlin] || c2>fHighThreshold[jlin]){
 593               fGoodInj[jpad][jlin]=1;
 594               break;
 595             }
 596           }
 597         }
 598       }
 599     }
 600   }
 601 }
 602 //______________________________________________________________________
 603 void AliITSOnlineSDDInjectors::FindCentroids(){
 604   // Computes the centroids (weighted mean) of teh injector pads
 605   for(Int_t jpad=0;jpad<kInjPads;jpad++){
 606     Int_t ian=GetAnodeNumber(jpad);
 607     for(Int_t jlin=0;jlin<kInjLines;jlin++){
 608       if(!fGoodInj[jpad][jlin]) continue;
 609       Double_t maxcont=0;
 610       Int_t ilmax=-1;
 611       for(Int_t jjj=fTbMin[jlin];jjj<fTbMax[jlin];jjj++){
 612         Double_t cont=fHisto->GetBinContent(jjj,ian+1);
 613         if(cont>maxcont){
 614           maxcont=cont;
 615           ilmax=jjj;
 616         }
 617       }
 618       Double_t intCont=0;
 619       Int_t jjj=ilmax;
 620       while(1){
 621         Double_t cont=fHisto->GetBinContent(jjj,ian+1);
 622         if(cont<fLowThreshold[jlin]) break;
 623         if(cont<fgkSaturation){
 624           fCentroid[jpad][jlin]+=cont*(Double_t)jjj;
 625           fRMSCentroid[jpad][jlin]+=cont*(Double_t)jjj*(Double_t)jjj;
 626           intCont+=cont;
 627         }
 628         jjj--;
 629       }
 630       jjj=ilmax+1;
 631       while(1){
 632         Double_t cont=fHisto->GetBinContent(jjj,ian+1);
 633         if(cont<fLowThreshold[jlin]) break;
 634         if(cont<fgkSaturation){
 635           fCentroid[jpad][jlin]+=cont*float(jjj);
 636           fRMSCentroid[jpad][jlin]+=cont*(Double_t)jjj*(Double_t)jjj;
 637           intCont+=cont;
 638         }
 639         jjj++;
 640       }
 641       if(intCont>0){
 642         fCentroid[jpad][jlin]/=intCont;
 643         fRMSCentroid[jpad][jlin]=TMath::Sqrt(fRMSCentroid[jpad][jlin]/intCont-fCentroid[jpad][jlin]*fCentroid[jpad][jlin])/TMath::Sqrt(intCont);
 644       }
 645       else{
 646         fCentroid[jpad][jlin]=0.;
 647         fRMSCentroid[jpad][jlin]=0.;
 648         fGoodInj[jpad][jlin]=0;
 649       }
 650       if(fRMSCentroid[jpad][jlin]==0) fGoodInj[jpad][jlin]=0;
 651     }
 652   }
 653 }
 654 //______________________________________________________________________
 655 void AliITSOnlineSDDInjectors::PrintInjectorStatus(){
 656   // Dumps the status bit of injector pads
 657   for(Int_t jpad=0;jpad<kInjPads;jpad++){
 658     printf("Line%d-Anode%d: %d %d %d\n",jpad,GetAnodeNumber(jpad),fGoodInj[jpad][0],fGoodInj[jpad][1],fGoodInj[jpad][2]);
 659   }
 660 }
 661 //______________________________________________________________________
 662 void AliITSOnlineSDDInjectors::PrintCentroids(){
 663   // Dumps the drift time centroids of injector pads
 664   for(Int_t jpad=0;jpad<kInjPads;jpad++){
 665     printf("Line%d-Anode%d: %f+-%f %f+-%f %f+-%f\n",jpad,GetAnodeNumber(jpad),fCentroid[jpad][0],fRMSCentroid[jpad][0],fCentroid[jpad][1],fRMSCentroid[jpad][1],fCentroid[jpad][2],fRMSCentroid[jpad][2]);
 666   }
 667 }
 668 //______________________________________________________________________
 669 void AliITSOnlineSDDInjectors::WriteToASCII(Int_t evNumb, UInt_t timeStamp, Int_t optAppend){
 670   // writes drift speed vs. anode fit parameters into an ASCII file
 671   // to be sent to FXS by the DA and processed by the SHUTTLE
 672
 673   TString outfilnam;
 674   outfilnam.Form("SDDinj_ddl%02dc%02d_sid%d.data",fDDL,fCarlos,fSide);
 675   FILE* outf;
 676   if(optAppend==0){
 677     outf=fopen(outfilnam.Data(),"w");
 678     fprintf(outf,"%d\n",fActualPolDegree);
 679   }
 680   else outf=fopen(outfilnam.Data(),"a");
 681   fprintf(outf,"%d   %d   ",evNumb,timeStamp);
 682   for(Int_t ic=0;ic<fPolDegree+1;ic++){
 683     fprintf(outf,"%G ",fParam[ic]);
 684   }
 685   fprintf(outf,"\n");
 686   fclose(outf);
 687 }
 688 //______________________________________________________________________
 689 TH1F* AliITSOnlineSDDInjectors::GetMeanDriftSpeedVsPadHisto() const{
 690   // Builds histogram of average drift speed vs. pad number
 691   TString hisnam;
 692   hisnam.Form("hdrsp%02dc%02ds%d",fDDL,fCarlos,fSide);
 693   TH1F* h=new TH1F(hisnam.Data(),"",kInjPads,-0.5,kInjPads-0.5);
 694   if(fNEvents>0){
 695     for(Int_t i=0;i<kInjPads;i++){
 696       h->SetBinContent(i+1,GetMeanDriftSpeed(i));
 697       Double_t rms=GetRMSDriftSpeed(i);
 698       Double_t err=0.;
 699       if(rms>0.) err=rms/TMath::Sqrt(fNEventsInPad[i]);
 700       h->SetBinError(i+1,err);
 701     }
 702   }
 703   return h;
 704 }
 705 //______________________________________________________________________
 706 Bool_t AliITSOnlineSDDInjectors::WriteToROOT(TFile *fil) const {
 707   // Writes the output histograms into a root file
 708   if(fil==0){
 709     AliWarning("Invalid pointer to ROOT file");
 710     return kFALSE;
 711   }
 712   TString hisnam;
 713   fil->cd();
 714   hisnam.Form("hdrsp%02dc%02ds%d",fDDL,fCarlos,fSide);
 715   TH1F hdsp(hisnam.Data(),"",kInjPads,-0.5,kInjPads-0.5);
 716   if(fNEvents==0){
 717     AliWarning("Zero analyzed events");
 718     return kFALSE;
 719   }
 720
 721   for(Int_t i=0;i<kInjPads;i++){
 722     hdsp.SetBinContent(i+1,GetMeanDriftSpeed(i));
 723     Double_t rms=GetRMSDriftSpeed(i);
 724     Double_t err=0.;
 725     if(rms>0.) err=rms/TMath::Sqrt(fNEventsInPad[i]);
 726     hdsp.SetBinError(i+1,err);
 727   }
 728   hdsp.Write();
 729   return kTRUE;
 730 }
 731 //______________________________________________________________________
 732 void AliITSOnlineSDDInjectors::WriteInjectorStatusToASCII(){
 733   // dump status of injectors encoded into UInt_t
 734   // 5 bits (value 0-31) to store number of pads with given status
 735   TString outfilnam;
 736   outfilnam.Form("SDDinj_ddl%02dc%02d_sid%d.data",fDDL,fCarlos,fSide);
 737   FILE* outf=fopen(outfilnam.Data(),"a");
 738   Int_t n[8]={0,0,0,0,0,0,0,0};
 739   for(Int_t jpad=fFirstPadForFit; jpad<=fLastPadForFit; jpad++){
 740     Int_t statusPad=GetInjPadStatus(jpad);
 741     ++n[statusPad];
 742   }
 743   UInt_t statusInj=0;
 744   statusInj+=(n[7]&0x1F)<<25; // bits 25-29: n. of pads with status 7
 745   statusInj+=(n[6]&0x1F)<<20; // bits 20-24: n. of pads with status 6
 746   statusInj+=(n[5]&0x1F)<<15; // bits 15-19: n. of pads with status 5
 747   statusInj+=(n[4]&0x1F)<<10; // bits 10-14: n. of pads with status 4
 748   statusInj+=(n[3]&0x1F)<<5;  // bits  5- 9: n. of pads with status 3
 749   statusInj+=(n[2]&0x1F);     // bits  0- 4: n. of pads with status 2
 750
 751   fprintf(outf,"-99 %u\n",statusInj); // -99 used in preprocessor to find line
 752                                       // with injector status info
 753   fclose(outf);
 754
 755 }