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