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