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