EMCAL/macros/testRRF.C

   1 /*
   2  Small macro for testing the raw response function used in
   3  AliEMCALRawUtils
   4
   5  J.L. Klay (Cal Poly)
   6
   7 */
   8 const Double_t fgTimeTrigger = 1.5E-06;
   9 const Double_t fTau = 2.35;
  10 const Double_t fOrder = 2.;
  11 const Double_t fgPedestal = 32;
  12 const Double_t fTimeBinWidth = 100E-9;
  13 const Int_t fgOverflow = 0x3FF;
  14 const Double_t fHighLowGain = 16.;
  15 const Double_t fRawFormatTimeMax = 256*fTimeBinWidth;
  16 const Double_t fgFEENoise = 3.;
  17
  18
  19 void testRRF(const Double_t damp = 20, const Double_t dtime = 1e-09) {
  20
  21   TH1I* adcHigh = new TH1I("adcHigh","adcHigh",256,0,255);
  22   TH1I* adcLow = new TH1I("adcLow","adcLow",256,0,255);
  23
  24   TF1* signalF = new TF1("signal",RawResponseFunction, 0, 256, 5);
  25   signalF->SetParameter(0,damp);
  26   signalF->SetParameter(1,(dtime+fgTimeTrigger)/fTimeBinWidth);
  27   signalF->SetParameter(2,fTau);
  28   signalF->SetParameter(3,fOrder);
  29   signalF->SetParameter(4,fgPedestal);
  30
  31   for(Int_t itime = 0; itime < 256; itime++) {
  32     Double_t signal = signalF->Eval(itime);
  33
  34     Double_t noise = gRandom->Gaus(0.,fgFEENoise);
  35     signal = sqrt(signal*signal + noise*noise);
  36
  37     cout << "itime = " << itime << " highgain = " << signal << " signalI = " << static_cast<Int_t>(signal +0.5);
  38
  39     if(static_cast<Int_t>(signal +0.5) > fgOverflow) {
  40       adcHigh->SetBinContent(itime+1,fgOverflow);
  41     } else {
  42       adcHigh->SetBinContent(itime+1,static_cast<Int_t>(signal +0.5));
  43     }
  44
  45     signal /= fHighLowGain;
  46
  47     cout << " lowgain = " << signal << " signalI = " << static_cast<Int_t>(signal +0.5) << endl;
  48
  49     if(static_cast<Int_t>(signal +0.5) > fgOverflow) {
  50       adcLow->SetBinContent(itime+1,fgOverflow);
  51     } else {
  52       adcLow->SetBinContent(itime+1,static_cast<Int_t>(signal +0.5));
  53     }
  54   }
  55
  56   TCanvas *c1 = new TCanvas("c1","c1",20,20,600,1000);
  57   c1->Divide(1,3);
  58   c1->cd(1);
  59   signalF->Draw();
  60   c1->cd(2);
  61   adcHigh->Draw();
  62   c1->cd(3);
  63   adcLow->Draw();
  64
  65 }
  66
  67 Double_t RawResponseFunction(Double_t *x, Double_t *par)
  68 {
  69   Double_t signal ;
  70   Double_t tau =par[2];
  71   Double_t N =par[3];
  72   Double_t ped = par[4];
  73   Double_t xx = ( x[0] - par[1] + tau ) / tau ;
  74
  75   if (xx <= 0)
  76     signal = ped ;
  77   else {
  78     signal = ped + par[0] * TMath::Power(xx , N) * TMath::Exp(N * (1 - xx )) ;
  79   }
  80   return signal ;
  81 }