fill the fraction of cell energy in cluster in case of unfolding ON
[u/mrichter/AliRoot.git] / 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 }