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