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
  23
  24 ClassImp(AliSignalProcesor)
  25
  26
  27 Double_t asymgauss(Double_t* x, Double_t* par)
  28 {
  29   // par[0] = normalization
  30   // par[1] = mean
  31   // par[2] = sigma
  32   // norm0  = 1
  33   // par[3] = lambda0
  34   // par[4] = norm1
  35   // par[5] = lambda1
  36   //
  37
  38   Double_t par1save = par[1];
  39   Double_t par2save = par[2];
  40   Double_t par3save = par[3];
  41   Double_t par5save = par[5];
  42   Double_t dx   = x[0]-par1save;
  43   //
  44   //
  45   Double_t  sigma2  = par2save*par2save;
  46   Double_t  sqrt2   = TMath::Sqrt(2.);
  47   if (-par3save*(dx-0.5*par3save*sigma2)>100) return 0;   // avoid overflow
  48   if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow
  49   if (TMath::Abs(par[4])>1) return 0;
  50   Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
  51     *(1-TMath::Erf((par3save*sigma2-dx)/(sqrt2*par2save)));
  52
  53   Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
  54     *(1-TMath::Erf((par5save*sigma2-dx)/(sqrt2*par2save)));
  55
  56
  57   return par[0]*(exp1+par[4]*exp2);
  58 }
  59
  60 Double_t asymgaussN(Double_t* x, Double_t* par)
  61 {
  62   // par[0] = normalization
  63   // par[1] = mean
  64   // par[2] = sigma
  65   // norm0  = 1
  66   // par[3] = lambda0
  67   // par[4] = norm1
  68   // par[5] = lambda1
  69   //
  70
  71   Double_t par1save = par[1];
  72   Double_t par2save = par[2];
  73   Double_t par3save = par[3];
  74   Double_t par5save = par[5];
  75   Double_t dx   = x[0]-par1save;
  76   //
  77   //
  78   Double_t  sigma2  = par2save*par2save;
  79   Double_t  sqrt2   = TMath::Sqrt(2.);
  80   if (-par3save*(dx-0.5*par3save*sigma2)>100) return 0;   // avoid overflow
  81   if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow
  82   if (TMath::Abs(par[4])>=1) return 0;
  83
  84   Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
  85     *0.5*(1-TMath::Erf((par3save*sigma2-dx)/(sqrt2*par2save)));
  86
  87   Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
  88     *0.5*(1-TMath::Erf((par5save*sigma2-dx)/(sqrt2*par2save)));
  89
  90
  91   return par[0]*(1.*exp1+par[4]*exp2)/(1.+par[4]);
  92 }
  93
  94
  95 TF1 * AliSignalProcesor::GetAsymGauss()
  96 {
  97   TF1 * f1 = new TF1("asymg",asymgaussN,-10,40,6);
  98   return f1;
  99 }
 100
 101
 102
 103 void AliSignalProcesor::SplineSmoother(Double_t *ampin, Double_t *ampout, Int_t n)
 104 {
 105   //
 106   //
 107   Float_t in[10000];
 108   Float_t out[10000];
 109   in[0] = ampin[0];
 110   in[1] = (ampin[0]+ampin[1])*0.5;
 111   in[2*(n-1)]    = ampin[n-1];
 112   in[2*(n-1)+1]  = ampin[n-1];
 113   //
 114   // add charge to the end
 115   for (Int_t i=0;i<10;i++){
 116     in[2*(n-1)+i]=ampin[n-1];
 117   }
 118
 119   //
 120   for (Int_t i=1;i<n-1;i++){
 121     in[2*i]    = ampin[i];
 122     in[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16.;
 123   }
 124   //
 125   out[0] = in[0];
 126   for (Int_t i=1;i<=2*n;i++){
 127     out[i]  = (9.*(in[i]+in[i+1])-(in[i-1]+in[i+2]))/16.;
 128   }
 129   //
 130   //
 131   for (int i=0;i<n;i++){
 132     ampout[i]      = out[2*i+1];
 133   }
 134 }
 135
 136
 137
 138
 139 void AliSignalProcesor::TailCancelationALTRO(Double_t *ampin, Double_t *ampout, Float_t K, Float_t L,
 140                       Int_t n)
 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(Double_t *ampin, Double_t *ampout, Float_t r, Float_t c,
 156                       Int_t n)
 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(Double_t *ampin, Double_t *ampout, Float_t lambda,
 171                       Int_t n)
 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(Double_t *ampin, Double_t *ampout, Float_t norm,
 209                                         Float_t lambda, Int_t n)
 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(Double_t *ampin, Double_t *ampout, Float_t lambda,
 251                       Int_t n)
 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 }