TPC/AliTPCRF1D.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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
  16 /*
  17 $Log$
  18 */
  19
  20 //-----------------------------------------------------------------------------
  21 //  $Header$
  22 //
  23 //
  24 //  Origin:   Marian Ivanov, Uni. of Bratislava, ivanov@fmph.uniba.sk
  25 //
  26 //  Declaration of class AliTPCRF1D
  27 //
  28 //-----------------------------------------------------------------------------
  29 #include "TMath.h"
  30 #include "AliTPCRF1D.h"
  31 #include "TF2.h"
  32 #include <iostream.h>
  33 #include "TCanvas.h"
  34 #include "TPad.h"
  35 #include "TStyle.h"
  36 #include "TH1.h"
  37
  38 extern TStyle * gStyle;
  39
  40 static Double_t funGauss(Double_t *x, Double_t * par)
  41 {
  42   return TMath::Exp(-(x[0]*x[0])/(2*par[0]*par[0]));
  43 }
  44
  45 static Double_t funCosh(Double_t *x, Double_t * par)
  46 {
  47   return 1/TMath::CosH(3.14159*x[0]/(2*par[0]));
  48 }
  49
  50 static Double_t funGati(Double_t *x, Double_t * par)
  51 {
  52   //par[1] = is equal to k3
  53   //par[0] is equal to pad wire distance
  54   Float_t K3=par[1];
  55   Float_t K3R=TMath::Sqrt(K3);
  56   Float_t K2=(TMath::Pi()/2)*(1-K3R/2.);
  57   Float_t K1=K2*K3R/(4*TMath::ATan(K3R));
  58   Float_t l=x[0]/par[0];
  59   Float_t tan2=TMath::TanH(K2*l);
  60   tan2*=tan2;
  61   Float_t res = K1*(1-tan2)/(1+K3*tan2);
  62   return res;
  63 }
  64
  65
  66
  67
  68 ///////////////////////////////////////////////////////////////////////////
  69 ///////////////////////////////////////////////////////////////////////////
  70 ///////////////////////////////////////////////////////////////////////////
  71 ///////////////////////////////////////////////////////////////////////////
  72
  73 AliTPCRF1D * gRF1D;
  74 ClassImp(AliTPCRF1D)
  75
  76
  77 AliTPCRF1D::AliTPCRF1D(Bool_t direct,Int_t np,Float_t step)
  78 {
  79   fDirect=direct;
  80   fNRF = np;
  81   fcharge = new Float_t[fNRF];
  82   fDSTEPM1=1./step;
  83   fSigma = 0;
  84   gRF1D = this;
  85   fGRF = 0;
  86   fkNorm = 0.5;
  87   fpadWidth = 3.5;
  88   forigsigma=0.;
  89   fOffset = 0.;
  90 }
  91
  92
  93 AliTPCRF1D::~AliTPCRF1D()
  94 {
  95   if (fcharge!=0) delete fcharge;
  96   if (fGRF !=0 ) fGRF->Delete();
  97 }
  98
  99 Float_t AliTPCRF1D::GetRF(Float_t xin)
 100 {
 101   //x xin DSTEP unit
 102   //return linear aproximation of RF
 103   Float_t x = TMath::Abs((xin-fOffset)*fDSTEPM1)+fNRF/2;
 104   Int_t i1=Int_t(x);
 105   if (x<0) i1-=1;
 106   Float_t res=0;
 107   if (i1+1<fNRF)
 108     res = fcharge[i1]*(Float_t(i1+1)-x)+fcharge[i1+1]*(x-Float_t(i1));
 109   return res;
 110 }
 111
 112 Float_t  AliTPCRF1D::GetGRF(Float_t xin)
 113 {
 114   if (fGRF != 0 )
 115     return fkNorm*fGRF->Eval(xin)/fInteg;
 116       else
 117     return 0.;
 118 }
 119
 120
 121 void AliTPCRF1D::SetParam( TF1 * GRF,Float_t padwidth,
 122                        Float_t kNorm, Float_t sigma)
 123 {
 124    fpadWidth = padwidth;
 125    fGRF = GRF;
 126    fkNorm = kNorm;
 127    if (sigma==0) sigma= fpadWidth/TMath::Sqrt(12.);
 128    forigsigma=sigma;
 129    fDSTEPM1 = 10/TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12);
 130    sprintf(fType,"User");
 131    //   Update();
 132 }
 133
 134
 135 void AliTPCRF1D::SetGauss(Float_t sigma, Float_t padWidth,
 136                       Float_t kNorm)
 137 {
 138   // char s[120];
 139   fpadWidth = padWidth;
 140   fkNorm = kNorm;
 141   if (fGRF !=0 ) fGRF->Delete();
 142   fGRF = new TF1("fun",funGauss,-5,5,2);
 143   funParam[0]=sigma;
 144   forigsigma=sigma;
 145   fGRF->SetParameters(funParam);
 146    fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12);
 147   //by default I set the step as one tenth of sigma
 148    //  Update();
 149   sprintf(fType,"Gauss");
 150 }
 151
 152 void AliTPCRF1D::SetCosh(Float_t sigma, Float_t padWidth,
 153                      Float_t kNorm)
 154 {
 155   //  char s[120];
 156   fpadWidth = padWidth;
 157   fkNorm = kNorm;
 158   if (fGRF !=0 ) fGRF->Delete();
 159   fGRF = new TF1("fun", funCosh, -5.,5.,2);
 160   funParam[0]=sigma;
 161   fGRF->SetParameters(funParam);
 162   forigsigma=sigma;
 163   fDSTEPM1 = 10./TMath::Sqrt(sigma*sigma+fpadWidth*fpadWidth/12);
 164   //by default I set the step as one tenth of sigma
 165   //  Update();
 166   sprintf(fType,"Cosh");
 167 }
 168
 169 void AliTPCRF1D::SetGati(Float_t K3, Float_t padDistance, Float_t padWidth,
 170                      Float_t kNorm)
 171 {
 172   //  char s[120];
 173   fpadWidth = padWidth;
 174   fkNorm = kNorm;
 175   if (fGRF !=0 ) fGRF->Delete();
 176   fGRF = new TF1("fun", funGati, -5.,5.,2);
 177   funParam[0]=padDistance;
 178   funParam[1]=K3;
 179   fGRF->SetParameters(funParam);
 180   forigsigma=padDistance;
 181   fDSTEPM1 = 10./TMath::Sqrt(padDistance*padDistance+fpadWidth*fpadWidth/12);
 182   //by default I set the step as one tenth of sigma
 183   //  Update();
 184   sprintf(fType,"Gati");
 185 }
 186
 187 void AliTPCRF1D::Draw(Float_t x1,Float_t x2,Int_t N)
 188 {
 189   char s[100];
 190   TCanvas  * c1 = new TCanvas("canRF","Pad response function",700,900);
 191   c1->cd();
 192   TPad * pad1 = new TPad("pad1RF","",0.05,0.55,0.95,0.95,21);
 193   pad1->Draw();
 194   TPad * pad2 = new TPad("pad2RF","",0.05,0.05,0.95,0.45,21);
 195   pad2->Draw();
 196
 197   sprintf(s,"RF response function for %1.2f cm pad width",
 198           fpadWidth);
 199   pad1->cd();
 200   TH1F * hRFo = new TH1F("hRFo","Original charge distribution",N+1,x1,x2);
 201   pad2->cd();
 202    gStyle->SetOptFit(1);
 203    gStyle->SetOptStat(0);
 204   TH1F * hRFc = new TH1F("hRFc",s,N+1,x1,x2);
 205   Float_t x=x1;
 206   Float_t y1;
 207   Float_t y2;
 208
 209   for (Float_t i = 0;i<N+1;i++)
 210     {
 211       x+=(x2-x1)/Float_t(N);
 212       y1 = GetRF(x);
 213       hRFc->Fill(x,y1);
 214       y2 = GetGRF(x);
 215       hRFo->Fill(x,y2);
 216     };
 217   pad1->cd();
 218   hRFo->Fit("gaus");
 219   pad2->cd();
 220   hRFc->Fit("gaus");
 221 }
 222
 223 void AliTPCRF1D::Update()
 224 {
 225   //initialize to 0
 226   for (Int_t i =0; i<fNRF;i++)  fcharge[i] = 0;
 227   if ( fGRF == 0 ) return;
 228   fInteg  = fGRF->Integral(-5*forigsigma,5*forigsigma,funParam,0.00001);
 229   if ( fInteg == 0 ) fInteg = 1;
 230   if (fDirect==kFALSE){
 231   //integrate charge over pad for different distance of pad
 232   for (Int_t i =0; i<fNRF;i++)
 233     {      //x in cm fpadWidth in cm
 234       Float_t x = (Float_t)(i-fNRF/2)/fDSTEPM1;
 235       Float_t x1=TMath::Max(x-fpadWidth/2,-5*forigsigma);
 236       Float_t x2=TMath::Min(x+fpadWidth/2,5*forigsigma);
 237       fcharge[i] =
 238         fkNorm*fGRF->Integral(x1,x2,funParam,0.0001)/fInteg;
 239     };
 240   }
 241   else{
 242     for (Int_t i =0; i<fNRF;i++)
 243       {      //x in cm fpadWidth in cm
 244         Float_t x = (Float_t)(i-fNRF/2)/fDSTEPM1;
 245         fcharge[i] = fkNorm*fGRF->Eval(x);
 246       };
 247   }
 248   fSigma = 0;
 249   Float_t sum =0;
 250   Float_t mean=0;
 251   for (Float_t  x =-fNRF/fDSTEPM1; x<fNRF/fDSTEPM1;x+=1/fDSTEPM1)
 252     {      //x in cm fpadWidth in cm
 253       Float_t weight = GetRF(x+fOffset);
 254       fSigma+=x*x*weight;
 255       mean+=x*weight;
 256       sum+=weight;
 257     };
 258   if (sum>0){
 259     mean/=sum;
 260     fSigma = TMath::Sqrt(fSigma/sum-mean*mean);
 261   }
 262   else fSigma=0;
 263 }
 264
 265 void AliTPCRF1D::Streamer(TBuffer &R__b)
 266 {
 267    // Stream an object of class AliTPC.
 268
 269    if (R__b.IsReading()) {
 270       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
 271       TObject::Streamer(R__b);
 272       //read pad parameters
 273       R__b >> fpadWidth;
 274       //read charge parameters
 275       R__b >> fType[0];
 276       R__b >> fType[1];
 277       R__b >> fType[2];
 278       R__b >> fType[3];
 279       R__b >> fType[4];
 280       R__b >> forigsigma;
 281       R__b >> fkNorm;
 282       R__b >> fK3X;
 283       R__b >> fPadDistance;
 284       R__b >> fInteg;
 285       R__b >> fOffset;
 286       //read functions
 287       if (fGRF!=0) {
 288         delete fGRF;
 289         fGRF=0;
 290       }
 291       if (strncmp(fType,"User",3)==0){
 292         fGRF= new TF1;
 293         R__b>>fGRF;
 294       }
 295
 296       if (strncmp(fType,"Gauss",3)==0)
 297         fGRF = new TF1("fun",funGauss,-5.,5.,4);
 298       if (strncmp(fType,"Cosh",3)==0)
 299         fGRF = new TF1("fun",funCosh,-5.,5.,4);
 300       if (strncmp(fType,"Gati",3)==0)
 301         fGRF = new TF1("fun",funGati,-5.,5.,4);
 302       R__b >>fDSTEPM1;
 303       R__b >>fNRF;
 304       R__b.ReadFastArray(fcharge,fNRF);
 305       R__b.ReadFastArray(funParam,5);
 306       if (fGRF!=0) fGRF->SetParameters(funParam);
 307
 308    } else {
 309       R__b.WriteVersion(AliTPCRF1D::IsA());
 310       TObject::Streamer(R__b);
 311       //write pad width
 312       R__b << fpadWidth;
 313       //write charge parameters
 314       R__b << fType[0];
 315       R__b << fType[1];
 316       R__b << fType[2];
 317       R__b << fType[3];
 318       R__b << fType[4];
 319       R__b << forigsigma;
 320       R__b << fkNorm;
 321       R__b << fK3X;
 322       R__b << fPadDistance;
 323       R__b << fInteg;
 324       R__b << fOffset;
 325       //write interpolation parameters
 326       if (strncmp(fType,"User",3)==0) R__b <<fGRF;
 327       R__b <<fDSTEPM1;
 328       R__b <<fNRF;
 329       R__b.WriteFastArray(fcharge,fNRF);
 330       R__b.WriteFastArray(funParam,5);
 331
 332
 333
 334    }
 335 }
 336