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