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 "AliITSOnlineSDDInjectors.h"
  16 #include <TH2F.h>
  17 #include <TGraphErrors.h>
  18 #include <TMath.h>
  19
  20 /* $Id$ */
  21
  22 ///////////////////////////////////////////////////////////////////
  23 //                                                               //
  24 // Implementation of the class used for SDD injector analysis    //
  25 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  26 //                                                               //
  27 ///////////////////////////////////////////////////////////////////
  28
  29 ClassImp(AliITSOnlineSDDInjectors)
  30
  31 const Float_t AliITSOnlineSDDInjectors::fgkSaturation = 1008.;
  32 const Float_t AliITSOnlineSDDInjectors::fgkJitterTB = 8.;
  33
  34 //______________________________________________________________________
  35   AliITSOnlineSDDInjectors::AliITSOnlineSDDInjectors():AliITSOnlineSDD(),fHisto(),fTbZero(0.),fParam(),fPolOrder(0),fMinDriftVel(0.),fMaxDriftVel(0.),fThreshold(0.)
  36 {
  37   // default constructor
  38   SetMinDriftVel();
  39   SetMaxDriftVel();
  40   SetRangeLine1();
  41   SetRangeLine2();
  42   SetRangeLine3();
  43   SetPositions();
  44   SetPolOrder();
  45   SetThreshold();
  46 }
  47 //______________________________________________________________________
  48 AliITSOnlineSDDInjectors::AliITSOnlineSDDInjectors(Int_t mod, Int_t sid):AliITSOnlineSDD(mod,sid),fHisto(),fTbZero(0.),fParam(),fPolOrder(0),fMinDriftVel(0.),fMaxDriftVel(0.),fThreshold(0.)
  49 {
  50 // standard constructor
  51   SetMinDriftVel();
  52   SetMaxDriftVel();
  53   SetRangeLine1();
  54   SetRangeLine2();
  55   SetRangeLine3();
  56   SetPositions();
  57   SetPolOrder();
  58   SetThreshold();
  59 }
  60 //______________________________________________________________________
  61 AliITSOnlineSDDInjectors::~AliITSOnlineSDDInjectors(){
  62   // Destructor
  63   if(fHisto) delete fHisto;
  64   if(fParam) delete [] fParam;
  65 }
  66 //______________________________________________________________________
  67 void AliITSOnlineSDDInjectors::SetPositions(){
  68   //
  69   Float_t kLinFromCenterUm[3]={31860.,17460.,660.};
  70   Float_t kAnodeFromCenterUm=35085;
  71   for(Int_t i=0;i<3;i++){
  72     fPosition[i]=kAnodeFromCenterUm-kLinFromCenterUm[i];
  73     fPosition[i]/=10000.; // from microns to cm
  74   }
  75 }
  76 //______________________________________________________________________
  77 void AliITSOnlineSDDInjectors::Reset(){
  78   //
  79   for(Int_t i=0;i<kNInjectors;i++){
  80     fDriftVel[i]=0.;
  81     fSigmaDriftVel[i]=0.;
  82   }
  83   for(Int_t i=0;i<kNInjectors;i++){
  84     for(Int_t j=0;j<3;j++){
  85       fGoodInj[i][j]=0;
  86       fCentroid[i][j]=0.;
  87       fRMSCentroid[i][j]=0.;
  88     }
  89   }
  90 }
  91 //______________________________________________________________________
  92 void AliITSOnlineSDDInjectors::AnalyzeEvent(TH2F* his){
  93   //
  94   Reset();
  95   fHisto=his;
  96   FindGoodInjectors();
  97   FindCentroids();
  98   CalcTimeBinZero();
  99   for(Int_t j=0;j<kNInjectors;j++) CalcDriftVelocity(j);
 100   FitDriftVelocityVsAnode();
 101 }
 102 //______________________________________________________________________
 103 TGraphErrors* AliITSOnlineSDDInjectors::GetLineGraph(Int_t jlin){
 104   //
 105   Float_t x[4],y[4],ex[4],ey[4];
 106   x[0]=0.;
 107   ex[0]=0.;
 108   y[0]=fTbZero;
 109   ey[0]=0.;
 110   for(Int_t i=0;i<3;i++){
 111     x[i+1]=fPosition[i];
 112     ex[i+1]=0.;
 113     y[i+1]=fCentroid[jlin][i];
 114     ey[i+1]=fRMSCentroid[jlin][i];
 115   }
 116   TGraphErrors *g=new TGraphErrors(4,x,y,ex,ey);
 117   return g;
 118 }
 119 //______________________________________________________________________
 120 Float_t AliITSOnlineSDDInjectors::GetDriftCoordinate(Float_t cAnode, Float_t cTimeBin){
 121   //
 122   Float_t vel=0;
 123   for(Int_t i=0;i<=fPolOrder;i++) vel+=fParam[i]*TMath::Power(cAnode,(Float_t)i);
 124   return vel*(cTimeBin-(fTbZero-fgkJitterTB))*25/1000.;
 125 }
 126 //______________________________________________________________________
 127 TGraphErrors* AliITSOnlineSDDInjectors::GetDriftVelocityGraph() const{
 128   //
 129   Int_t ipt=0;
 130   TGraphErrors *g=new TGraphErrors(0);
 131   for(Int_t i=0;i<kNInjectors;i++){
 132     if(fDriftVel[i]>0){
 133       g->SetPoint(ipt,GetAnodeNumber(i),fDriftVel[i]);
 134       g->SetPointError(ipt,0,fSigmaDriftVel[i]);
 135       ipt++;
 136     }
 137   }
 138   return g;
 139 }
 140 //______________________________________________________________________
 141 void AliITSOnlineSDDInjectors::CalcTimeBinZero(){
 142   //
 143   Float_t tzero=0.,intCont=0.;
 144   for(Int_t ian=0;ian<fgkNAnodes;ian++){
 145     for(Int_t itb=1;itb<fTbMin[0];itb++){
 146       Float_t cont=fHisto->GetBinContent(itb,ian+1);
 147       if(cont>fThreshold){
 148         tzero+=cont*float(itb);
 149         intCont+=cont;
 150       }
 151     }
 152   }
 153   if(intCont>0) fTbZero=tzero/intCont;
 154 }
 155 //______________________________________________________________________
 156 void AliITSOnlineSDDInjectors::FitDriftVelocityVsAnode(){
 157   const Int_t nn=fPolOrder+1;
 158   Float_t **mat = new Float_t*[nn];
 159   for(Int_t i=0; i < nn; i++) mat[i] = new Float_t[nn];
 160   Float_t *vect = new Float_t[nn];
 161   for(Int_t k1=0;k1<nn;k1++){
 162     vect[k1]=0;
 163     for(Int_t k2=0;k2<nn;k2++){
 164       mat[k1][k2]=0;
 165       for(Int_t n=0; n<kNInjectors;n++){
 166         Float_t x=(Float_t)GetAnodeNumber(n);
 167         if(fDriftVel[n]>0) mat[k1][k2]+=TMath::Power(x,k1+k2)/TMath::Power(fSigmaDriftVel[n],2);
 168       }
 169     }
 170   }
 171   for(Int_t k1=0;k1<nn;k1++){
 172     for(Int_t n=0; n<kNInjectors;n++){
 173       Float_t x=(Float_t)GetAnodeNumber(n);
 174       if(fDriftVel[n]>0) vect[k1]+=fDriftVel[n]*TMath::Power(x,k1)/TMath::Power(fSigmaDriftVel[n],2);
 175     }
 176   }
 177   Int_t *iPivot = new Int_t[nn];
 178   Int_t *indxR = new Int_t[nn];
 179   Int_t *indxC = new Int_t[nn];
 180   for(Int_t i=0;i<nn;i++) iPivot[i]=0;
 181   Int_t iCol=-1,iRow=-1;
 182   for(Int_t i=0;i<nn;i++){
 183     Float_t big=0.;
 184     for(Int_t j=0;j<nn;j++){
 185       if(iPivot[j]!=1){
 186         for(Int_t k=0;k<nn;k++){
 187            if(iPivot[k]==0){
 188              if(TMath::Abs(mat[j][k])>=big){
 189                big=TMath::Abs(mat[j][k]);
 190                iRow=j;
 191                iCol=k;
 192              }
 193           }
 194         }
 195       }
 196     }
 197     iPivot[iCol]++;
 198     Float_t aux;
 199     if(iRow!=iCol){
 200       for(Int_t l=0;l<nn;l++){
 201         aux=mat[iRow][l];
 202         mat[iRow][l]=mat[iCol][l];
 203         mat[iCol][l]=aux;
 204       }
 205       aux=vect[iRow];
 206       vect[iRow]=vect[iCol];
 207       vect[iCol]=aux;
 208     }
 209     indxR[i]=iRow;
 210     indxC[i]=iCol;
 211     if(mat[iCol][iCol]==0) break;
 212     Float_t pivinv=1./mat[iCol][iCol];
 213     mat[iCol][iCol]=1;
 214     for(Int_t l=0;l<nn;l++) mat[iCol][l]*=pivinv;
 215     vect[iCol]*=pivinv;
 216     for(Int_t m=0;m<nn;m++){
 217       if(m!=iCol){
 218         aux=mat[m][iCol];
 219         mat[m][iCol]=0;
 220         for(Int_t n=0;n<nn;n++) mat[m][n]-=mat[iCol][n]*aux;
 221         vect[m]-=vect[iCol]*aux;
 222       }
 223     }
 224   }
 225   delete [] iPivot;
 226   delete [] indxR;
 227   delete [] indxC;
 228
 229   if(fParam) delete [] fParam;
 230   fParam=new Float_t[nn];
 231   for(Int_t i=0; i<nn;i++)fParam[i]=vect[i];
 232
 233   for(Int_t i=0; i < nn; i++) delete [] mat[i];
 234   delete [] mat;
 235   delete [] vect;
 236 }
 237 //______________________________________________________________________
 238 void AliITSOnlineSDDInjectors::CalcDriftVelocity(Int_t jlin){
 239   //
 240   Float_t sumY=0,sumX=0,sumXX=0,sumYY=0.,sumXY=0,sumWEI=0.;
 241   Int_t npt=0;
 242   Float_t y[3],ey[3];
 243   Float_t tzero=0,erry=0;
 244   for(Int_t i=0;i<3;i++){
 245     y[i]=fCentroid[jlin][i];
 246     ey[i]=fRMSCentroid[jlin][i];
 247   }
 248   for(Int_t i=0;i<3;i++){
 249     if(fGoodInj[jlin][i]){
 250       sumY+=y[i]/ey[i]/ey[i];
 251       sumX+=fPosition[i]/ey[i]/ey[i];
 252       sumXX+=fPosition[i]*fPosition[i]/ey[i]/ey[i];
 253       sumYY+=y[i]*y[i]/ey[i]/ey[i];
 254       sumXY+=fPosition[i]*y[i]/ey[i]/ey[i];
 255       sumWEI+=1./ey[i]/ey[i];
 256       tzero=fTbZero/ey[i]/ey[i];
 257       erry=ey[i];
 258       npt++;
 259     }
 260   }
 261   Float_t vel=0,evel=0;
 262   if(npt>1){
 263     Float_t slope=(sumWEI*sumXY-sumY*sumX)/(sumWEI*sumXX-sumX*sumX);
 264     Float_t eslope=TMath::Sqrt(sumWEI/(sumWEI*sumXX-sumX*sumX));
 265     vel=1./slope*10000./25.;// micron/ns
 266     evel=eslope/slope/slope*10000./25.;// micron/ns
 267   }
 268   if(npt==1){
 269     Float_t slope=(sumY-tzero)/sumX;
 270     Float_t eslope=erry/sumX;
 271     vel=1./slope*10000./25.;// micron/ns
 272     evel=eslope/slope/slope*10000./25.;// micron/ns
 273   }
 274   if(vel>fMaxDriftVel||vel<fMinDriftVel){
 275     vel=0.;
 276     evel=0.;
 277   }
 278   fDriftVel[jlin]=vel;
 279   fSigmaDriftVel[jlin]=evel;
 280 }
 281 //______________________________________________________________________
 282 Int_t AliITSOnlineSDDInjectors::GetAnodeNumber(Int_t iInjLine) const{
 283   //
 284   Int_t ian=-1;
 285   if(iInjLine>32) return ian;
 286   if(!fSide){
 287     ian=iInjLine*8;
 288     if(iInjLine==32) ian--;
 289   }else{
 290     ian=iInjLine*8-1;
 291     if(iInjLine==0) ian=0;
 292   }
 293   return ian;
 294 }
 295
 296 //______________________________________________________________________
 297 void AliITSOnlineSDDInjectors::FindGoodInjectors(){
 298   //
 299   for(Int_t iii=0;iii<kNInjectors;iii++){
 300     Int_t ian=GetAnodeNumber(iii);
 301     for(Int_t ninj=0;ninj<3;ninj++){
 302       for(Int_t jjj=fTbMin[ninj];jjj<fTbMax[ninj];jjj++){
 303         Float_t c1=fHisto->GetBinContent(jjj,ian+1);
 304         Float_t c2=fHisto->GetBinContent(jjj+1,ian+1);
 305         Float_t c3=fHisto->GetBinContent(jjj+2,ian+1);
 306         if(c1>fThreshold && c2>fThreshold && c3>fThreshold){
 307           fGoodInj[iii][ninj]=1;
 308           break;
 309         }
 310       }
 311       //      for(Int_t jjj=fTbMin[ninj];jjj<fTbMax[ninj];jjj++){
 312       //        Float_t c1=fHisto->GetBinContent(jjj,ian+1);
 313       //        if(c1>=fgkSaturation){
 314       //          fGoodInj[iii][ninj]=0;
 315       //          break;
 316       //        }
 317       //      }
 318     }
 319   }
 320 }
 321 //______________________________________________________________________
 322 void AliITSOnlineSDDInjectors::FindCentroids(){
 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       if(!fGoodInj[iii][ninj]) continue;
 328       Float_t maxcont=0;
 329       Int_t ilmax=-1;
 330       for(Int_t jjj=fTbMin[ninj];jjj<fTbMax[ninj];jjj++){
 331         Float_t cont=fHisto->GetBinContent(jjj,ian+1);
 332         if(cont>maxcont){
 333           maxcont=cont;
 334           ilmax=jjj;
 335         }
 336       }
 337       Float_t intCont=0;
 338       Int_t jjj=ilmax;
 339       while(1){
 340         Float_t cont=fHisto->GetBinContent(jjj,ian+1);
 341         if(cont<fThreshold) break;
 342         if(cont<fgkSaturation){
 343           fCentroid[iii][ninj]+=cont*(Float_t)jjj;
 344           fRMSCentroid[iii][ninj]+=cont*TMath::Power((Float_t)jjj,2);
 345           intCont+=cont;
 346         }
 347         jjj--;
 348       }
 349       jjj=ilmax+1;
 350       while(1){
 351         Float_t cont=fHisto->GetBinContent(jjj,ian+1);
 352         if(cont<fThreshold) break;
 353         if(cont<fgkSaturation){
 354           fCentroid[iii][ninj]+=cont*float(jjj);
 355           fRMSCentroid[iii][ninj]+=cont*TMath::Power((Float_t)jjj,2);
 356           intCont+=cont;
 357         }
 358         jjj++;
 359       }
 360       if(intCont>0){
 361         fCentroid[iii][ninj]/=intCont;
 362         fRMSCentroid[iii][ninj]=TMath::Sqrt(fRMSCentroid[iii][ninj]/intCont-fCentroid[iii][ninj]*fCentroid[iii][ninj]);
 363       }
 364       else{
 365         fCentroid[iii][ninj]=0.;
 366         fRMSCentroid[iii][ninj]=0.;
 367         fGoodInj[iii][ninj]=0;
 368       }
 369       if(fRMSCentroid[iii][ninj]==0) fGoodInj[iii][ninj]=0;
 370     }
 371   }
 372 }
 373 //______________________________________________________________________
 374 void AliITSOnlineSDDInjectors::PrintInjMap(){
 375   //
 376   for(Int_t iii=0;iii<kNInjectors;iii++){
 377     printf("Line%d-Anode%d: %d %d %d\n",iii,GetAnodeNumber(iii),fGoodInj[iii][0],fGoodInj[iii][1],fGoodInj[iii][2]);
 378   }
 379 }
 380 //______________________________________________________________________
 381 void AliITSOnlineSDDInjectors::PrintCentroids(){
 382   //
 383   for(Int_t iii=0;iii<kNInjectors;iii++){
 384     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]);
 385   }
 386 }
 387 //______________________________________________________________________
 388 void AliITSOnlineSDDInjectors::WriteToASCII(Int_t evNumb, UInt_t timeStamp, Int_t optAppend){
 389   //
 390   Char_t outfilnam[100];
 391   sprintf(outfilnam,"SDDinj_mod%03d_sid%d.data",fModuleId,fSide);
 392   FILE* outf;
 393   if(optAppend==0) outf=fopen(outfilnam,"w");
 394   else outf=fopen(outfilnam,"a");
 395   fprintf(outf,"%d   %d   ",evNumb,timeStamp);
 396   for(Int_t ic=0;ic<fPolOrder+1;ic++){
 397     fprintf(outf,"%G ",fParam[ic]);
 398   }
 399   fprintf(outf,"\n");
 400   fclose(outf);
 401 }