TPC/AliTPCRF1D.cxx

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