STEER/AliSignalProcesor.cxx

   1 #include "TF1.h"
   2 #include "AliSignalProcesor.h"
   3
   4
   5 ClassImp(AliSignalProcesor)
   6
   7
   8 Double_t asymgauss(Double_t* x, Double_t* par)
   9 {
  10   // par[0] = normalization
  11   // par[1] = mean
  12   // par[2] = sigma
  13   // norm0  = 1
  14   // par[3] = lambda0
  15   // par[4] = norm1
  16   // par[5] = lambda1
  17   //
  18
  19   Double_t par1save = par[1];
  20   Double_t par2save = par[2];
  21   Double_t par3save = par[3];
  22   Double_t par5save = par[5];
  23   Double_t dx   = x[0]-par1save;
  24   //
  25   //
  26   Double_t  sigma2  = par2save*par2save;
  27   Double_t  sqrt2   = TMath::Sqrt(2.);
  28   if (-par3save*(dx-0.5*par3save*sigma2)>100) return 0;   // avoid overflow
  29   if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow
  30   if (TMath::Abs(par[4])>1) return 0;
  31   Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
  32     *(1-TMath::Erf((par3save*sigma2-dx)/(sqrt2*par2save)));
  33
  34   Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
  35     *(1-TMath::Erf((par5save*sigma2-dx)/(sqrt2*par2save)));
  36
  37
  38   return par[0]*(exp1+par[4]*exp2);
  39 }
  40
  41 Double_t asymgaussN(Double_t* x, Double_t* par)
  42 {
  43   // par[0] = normalization
  44   // par[1] = mean
  45   // par[2] = sigma
  46   // norm0  = 1
  47   // par[3] = lambda0
  48   // par[4] = norm1
  49   // par[5] = lambda1
  50   //
  51
  52   Double_t par1save = par[1];
  53   Double_t par2save = par[2];
  54   Double_t par3save = par[3];
  55   Double_t par5save = par[5];
  56   Double_t dx   = x[0]-par1save;
  57   //
  58   //
  59   Double_t  sigma2  = par2save*par2save;
  60   Double_t  sqrt2   = TMath::Sqrt(2.);
  61   if (-par3save*(dx-0.5*par3save*sigma2)>100) return 0;   // avoid overflow
  62   if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow
  63   if (TMath::Abs(par[4])>=1) return 0;
  64
  65   Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
  66     *0.5*(1-TMath::Erf((par3save*sigma2-dx)/(sqrt2*par2save)));
  67
  68   Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
  69     *0.5*(1-TMath::Erf((par5save*sigma2-dx)/(sqrt2*par2save)));
  70
  71
  72   return par[0]*(1.*exp1+par[4]*exp2)/(1.+par[4]);
  73 }
  74
  75
  76 TF1 * AliSignalProcesor::GetAsymGauss()
  77 {
  78   TF1 * f1 = new TF1("asymg",asymgaussN,-10,40,6);
  79   return f1;
  80 }
  81
  82
  83
  84 void AliSignalProcesor::SplineSmoother(Double_t *ampin, Double_t *ampout, Int_t n)
  85 {
  86   //
  87   //
  88   Float_t in[10000];
  89   Float_t out[10000];
  90   in[0] = ampin[0];
  91   in[1] = (ampin[0]+ampin[1])*0.5;
  92   in[2*(n-1)]    = ampin[n-1];
  93   in[2*(n-1)+1]  = ampin[n-1];
  94   //
  95   // add charge to the end
  96   for (Int_t i=0;i<10;i++){
  97     in[2*(n-1)+i]=ampin[n-1];
  98   }
  99
 100   //
 101   for (Int_t i=1;i<n-1;i++){
 102     in[2*i]    = ampin[i];
 103     in[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16.;
 104   }
 105   //
 106   out[0] = in[0];
 107   for (Int_t i=1;i<=2*n;i++){
 108     out[i]  = (9.*(in[i]+in[i+1])-(in[i-1]+in[i+2]))/16.;
 109   }
 110   //
 111   //
 112   for (int i=0;i<n;i++){
 113     ampout[i]      = out[2*i+1];
 114   }
 115 }
 116
 117
 118
 119
 120 void AliSignalProcesor::TailCancelationALTRO(Double_t *ampin, Double_t *ampout, Float_t K, Float_t L,
 121                       Int_t n)
 122 {
 123   //
 124   // ALTRO
 125   Float_t temp;
 126   ampout[0]  = ampin[0];
 127   temp = ampin[0];
 128   for (int i=1;i<n;i++){
 129     ampout[i]   = ampin[i]   + (K-L)*temp;
 130     temp        = ampin[i]   +  K*temp;
 131   }
 132 }
 133
 134 //
 135 //
 136 void AliSignalProcesor::TailCancelationTRD(Double_t *ampin, Double_t *ampout, Float_t r, Float_t c,
 137                       Int_t n)
 138 {
 139   //TRD
 140   //
 141   Double_t reminder=0;
 142   //
 143   for (Int_t i=0; i<n; i++){
 144     ampout[i] = ampin[i]-reminder;
 145     //
 146     reminder  = r*(reminder+c*ampout[i]);
 147   }
 148
 149 }
 150
 151 void AliSignalProcesor::TailMaker(Double_t *ampin, Double_t *ampout, Float_t lambda,
 152                       Int_t n)
 153 {
 154
 155   Double_t l = TMath::Exp(-lambda);
 156   //
 157   Float_t temp=0;
 158   for (Int_t i=n-1; i>0; i--){
 159     ampout[i] = ampin[i]+temp;
 160     //
 161     temp  = l*(temp+ampin[i]);
 162   }
 163 }
 164
 165 void AliSignalProcesor::TailCancelationALTRO1(Double_t *ampin, Double_t *ampout, Float_t norm,
 166                                             Float_t lambda, Int_t n)
 167 {
 168
 169   Double_t l = TMath::Exp(-lambda);
 170   Double_t k = l*(1.-norm*lambda);
 171
 172   return TailCancelationALTRO(ampin,ampout,k,l,n);
 173 }
 174
 175
 176 void AliSignalProcesor::TailCancelationTRD1(Double_t *ampin, Double_t *ampout, Float_t norm,
 177                                           Float_t lambda, Int_t n)
 178 {
 179   //
 180   //
 181   Double_t r = TMath::Exp(-lambda);
 182   Double_t c = norm*lambda;
 183   return TailCancelationTRD(ampin,ampout,r,c,n);
 184 }
 185
 186
 187
 188
 189 void AliSignalProcesor::TailCancelationMI(Double_t *ampin, Double_t *ampout, Float_t norm,
 190                                         Float_t lambda, Int_t n)
 191 {
 192
 193   Double_t L = TMath::Exp(-lambda*0.5);
 194   Double_t K = L*(1.-norm*lambda*0.5);
 195   //
 196   //
 197   Float_t in[10000];
 198   Float_t out[10000];
 199   for (Int_t i=0;i<n*2+20;i++){
 200     in[i] = 0;
 201     out[i]= 0;
 202   }
 203   in[0] = ampin[0];
 204   in[1] = (ampin[0]+ampin[1])*0.5;
 205   in[2*(n-1)]    = ampin[n-1];
 206   in[2*(n-1)+1]  = ampin[n-1];
 207   //
 208   for (Int_t i=1;i<n-2;i++){
 209     in[2*i]    = ampin[i];
 210     in[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16;
 211   }
 212   //
 213   Float_t temp;
 214   out[0]     = in[0];
 215   temp       = in[0];
 216   for (int i=1;i<=2*n;i++){
 217     out[i]      = in[i]   + (K-L)*temp;
 218     temp        = in[i]   +  K*temp;
 219   }
 220   //
 221   //
 222   for (int i=0;i<n;i++){
 223     ampout[i]      = out[2*i+1];
 224   }
 225 }
 226
 227
 228
 229
 230
 231 void AliSignalProcesor::TailMakerSpline(Double_t *ampin, Double_t *ampout, Float_t lambda,
 232                       Int_t n)
 233 {
 234
 235   Double_t l = TMath::Exp(-lambda*0.5);
 236   //
 237   //
 238   Float_t in[10000];
 239   Float_t out[10000];
 240   for (Int_t i=0;i<n*2+20;i++){
 241     in[i] = 0;
 242     out[i]= 0;
 243   }
 244   in[0] = ampin[0];
 245   in[1] = (ampin[0]+ampin[1])*0.5;
 246   in[2*(n-1)]    = ampin[n-1];
 247   in[2*(n-1)+1]  = ampin[n-1];
 248   //
 249   // add charge to the end
 250   for (Int_t i=0;i<10;i++){
 251     in[2*(n-1)+i]=ampin[n-1];
 252   }
 253
 254   //
 255   for (Int_t i=1;i<n-2;i++){
 256     in[2*i]    = ampin[i];
 257     in[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16;
 258   }
 259   //
 260   //
 261   Float_t temp;
 262   out[2*n]     = in[2*n];
 263   temp         = 0;
 264   for (int i=2*n;i>=0;i--){
 265     out[i]      = in[i]   + temp;
 266     temp        = l*(temp+in[i]);
 267   }
 268   //
 269   //
 270   for (int i=0;i<n;i++){
 271     ampout[i]      = out[2*i+1];
 272   }
 273 }