For rare clusters with equal charge pads do just CoG().
[u/mrichter/AliRoot.git] / T0 / AliT0CalibWalk.cxx
CommitLineData
a2ad8166 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///////////////////////////////////////////////////////////////////////////////
19// //
20// class T0 walk correction Alla Maevskaya alla@inr.ru 20.11.2007 //
21// //
22///////////////////////////////////////////////////////////////////////////////
23
24#include "AliT0CalibWalk.h"
a2ad8166 25#include "AliLog.h"
a2ad8166 26
27#include <TObjArray.h>
28#include <TGraph.h>
29#include <TFile.h>
a2ad8166 30#include <TH2F.h>
31#include <TMath.h>
32#include <TSystem.h>
33#include <Riostream.h>
34#include <TSpectrum.h>
a2ad8166 35#include <TProfile.h>
36
a2ad8166 37
38ClassImp(AliT0CalibWalk)
39
40//________________________________________________________________
41 AliT0CalibWalk::AliT0CalibWalk(): TNamed(),
42 fWalk(0),
43 fAmpLEDRec(0)
44{
45 //
46}
47
48//________________________________________________________________
49AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
50 fWalk(0),
51 fAmpLEDRec(0)
52{
53 TString namst = "Calib_";
54 namst += name;
55 SetName(namst.Data());
56 SetTitle(namst.Data());
57
58}
59
60//________________________________________________________________
61AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
62 TNamed(calibda),
63 fWalk(0),
64 fAmpLEDRec(0)
65
66{
67// copy constructor
68 SetName(calibda.GetName());
69 SetTitle(calibda.GetName());
70
71
72}
73
74//________________________________________________________________
75AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
76{
77// assignment operator
78 SetName(calibda.GetName());
79 SetTitle(calibda.GetName());
80
81 return *this;
82}
83
84//________________________________________________________________
85AliT0CalibWalk::~AliT0CalibWalk()
86{
87 //
88}
89//________________________________________________________________
90void AliT0CalibWalk::SetWalk(Int_t ipmt)
91{
7576945f 92 //read QTC walk graph from external file
a2ad8166 93
94 Int_t mv, ps;
95 Int_t x[70000], y[70000], index[70000];
96 Float_t time[10000],amplitude[10000];
97 string buffer;
98 Bool_t down=false;
99
100 const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-Amp.txt");
101 ifstream inFile(filename);
102 if(!inFile) AliError(Form("Cannot open file %s !",filename));
103
104 Int_t i=0;
105 while(getline(inFile,buffer)){
106 inFile >> ps >> mv;
107
108 x[i]=ps; y[i]=mv;
109 i++;
110 }
111 inFile.close();
112 cout<<" number of data "<<i<<endl;
113
114 TMath::Sort(i, y, index,down);
115 Int_t amp=0, iin=0, isum=0, sum=0;
116 Int_t ind=0;
117 for (Int_t ii=0; ii<i; ii++)
118 {
119 ind=index[ii];
120 if(y[ind] == amp)
121 {
122 sum +=x[ind];
123 iin++;
124 }
125 else
126 {
127 if(iin>0)
128 time[isum] = Float_t (sum/(iin));
129 else
130 time[isum] =Float_t (x[ind]);
131 amplitude[isum] = Float_t (amp);
132 amp=y[ind];
133 iin=0;
134 isum++;
135 sum=0;
136 }
137 }
138
139 inFile.close();
140
141 TGraph* gr = new TGraph(isum, amplitude, time);
142 fWalk.AddAtAndExpand(gr,ipmt);
143}
144
145//________________________________________________________________
146
147void AliT0CalibWalk::SetAmpLEDRec(Int_t ipmt)
148{
7576945f 149 // read LED walk from external file
a2ad8166 150 Float_t mv, ps;
151 Float_t x[100], y[100];
152 string buffer;
153
154 const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-LED.txt");
155 ifstream inFile(filename);
156 if(!inFile) {AliError(Form("Cannot open file %s !",filename));}
157
158 inFile >> mv>>ps;
159 Int_t i=0;
160
161 while(getline(inFile,buffer)){
162 x[i]=mv; y[i]=ps;
163 inFile >> mv >> ps;
164 i++;
165 }
166 inFile.close();
167 Float_t y1[100], x1[100];
168 for (Int_t ir=0; ir<i; ir++){
169 y1[ir]=y[i-ir]; x1[ir]=x[i-ir];}
170 TGraph* gr = new TGraph(i,y1,x1);
171 fAmpLEDRec.AddAtAndExpand(gr,ipmt);
172
173}
174
175//________________________________________________________________
176
177void AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
178{
7576945f 179 //make walk corerction for preprocessor
180
a2ad8166 181 TFile *gFile = TFile::Open(laserFile);
700e6b36 182
a2ad8166 183 Int_t npeaks = 20;
083f0bec 184 Double_t sigma=3.;
700e6b36 185 Bool_t down = false;
2e6a5ee0 186
a2ad8166 187 Int_t index[20];
2e6a5ee0 188 Char_t buf1[20], buf2[20],buf3[20];
700e6b36 189
190
191 for (Int_t i=0; i<24; i++)
a2ad8166 192 {
700e6b36 193 if(i>11) { sprintf(buf1,"T0_A_%i_CFD",i+1-12); }
6b7d32f7 194
700e6b36 195 if(i<12) { sprintf(buf1,"T0_C_%i_CFD",i+1); }
f35a2c6c 196
2e6a5ee0 197 sprintf(buf2,"CFDvsQTC%i",i+1);
198 sprintf(buf3,"CFDvsLED%i",i+1);
199 // cout<<buf1<<" "<<buf2<<" "<<buf3<<endl;
083f0bec 200 TH2F *qtcVScfd = (TH2F*) gFile->Get(buf2);
2e6a5ee0 201 // qtcVScfd->Print();
083f0bec 202 TH2F *ledVScfd = (TH2F*) gFile->Get(buf3);
700e6b36 203 TH1F *cfd = (TH1F*) gFile->Get(buf1);
2e6a5ee0 204
700e6b36 205 TSpectrum *s = new TSpectrum(2*npeaks,1);
2e6a5ee0 206 Int_t nfound = s->Search(cfd,sigma,"goff",0.1);
207 // cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
208 if(nfound!=0){
700e6b36 209 Float_t *xpeak = s->GetPositionX();
210 TMath::Sort(nfound, xpeak, index,down);
211 Float_t xp = xpeak[index[0]];
212 Int_t xbin = cfd->GetXaxis()->FindBin(xp);
213 Float_t yp = cfd->GetBinContent(xbin);
2e6a5ee0 214 // cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
215 Float_t hmax = xp+20*sigma;
216 Float_t hmin = xp-20*sigma;
217 }
083f0bec 218 //QTC
219 Int_t nbins= qtcVScfd->GetXaxis()->GetNbins();
083f0bec 220 TProfile *prY = qtcVScfd->ProfileX();
2e6a5ee0 221 Float_t hmin = qtcVScfd->GetYaxis()->GetXmin();
222 Float_t hmax = qtcVScfd->GetYaxis()->GetXmax();
083f0bec 223 prY->SetMaximum(hmax);
224 prY->SetMinimum(hmin);
2e6a5ee0 225 Int_t np=0; //=nbins/5;
226 Double_t *xx = new Double_t[nbins];
227 Double_t *yy = new Double_t[nbins];
700e6b36 228 Int_t ng=0;
229 Double_t yg=0;
230 for (Int_t ip=1; ip<nbins; ip++)
231 {
2e6a5ee0 232
233 if(ip%5 != 0 )
234 {
235 // if (prY->GetBinContent(ip) !=0)
236 if (prY->GetBinContent(ip) >hmin && prY->GetBinContent(ip)<hmax){
083f0bec 237 yg +=prY->GetBinContent(ip);
2e6a5ee0 238 ng++;
239 // cout<<" ip "<<ip<<" Y "<<yg<<endl;
240 }
241 }
700e6b36 242 else {
2e6a5ee0 243 if (/*prY->GetBinContent(ip) >hmin && prY->GetBinContent(ip)<hmax &&*/ ng>0){
244 // xx[ip/5] = Float_t (prY->GetBinCenter(ip));
245 // yy[ip/5] = yg/ng;
246 xx[np] = Float_t (prY->GetBinCenter(ip));
247 yy[np] = yg/ng;
248 yg=0;
249 ng=0;
250 // cout<<" ip "<<ip<<" np "<<np<<" XX "<<xx[np]<<" YY "<<yy[np]<<endl;
251 np++;
252 }
700e6b36 253 }
254 }
2e6a5ee0 255
700e6b36 256 TGraph *gr = new TGraph(np,xx,yy);
257 gr->SetMinimum(hmin);
258 gr->SetMaximum(hmax);
259 fWalk.AddAtAndExpand(gr,i);
700e6b36 260
083f0bec 261 //LED
262 Int_t nbinsled= ledVScfd->GetXaxis()->GetNbins();
700e6b36 263 cout<<" nbins led "<<nbinsled<<endl;
083f0bec 264 TProfile *prledY = ledVScfd->ProfileX();
700e6b36 265
2e6a5ee0 266 Int_t npled = 0; //nbinsled/5;
267 Double_t *xxled = new Double_t[nbinsled];
268 Double_t *yyled = new Double_t[nbinsled];
700e6b36 269 Int_t ngled=0;
270 Double_t ygled=0;
271 for (Int_t ip=1; ip<nbinsled; ip++)
272 {
2e6a5ee0 273 if(ip%5 != 0 ) {
083f0bec 274 if (prledY->GetBinContent(ip) !=0)
275 ygled +=prledY->GetBinContent(ip);
700e6b36 276 ngled++;}
277 else {
2e6a5ee0 278 xxled[npled] = Float_t (prledY->GetBinCenter(ip));
279 yyled[npled] = ygled/ngled;
700e6b36 280 ygled=0;
281 ngled=0;
2e6a5ee0 282 npled++;
700e6b36 283 }
284 }
2e6a5ee0 285
286 TGraph *grled = new TGraph(npled,xxled,yyled);
700e6b36 287 grled->SetMinimum(hmin);
288 grled->SetMaximum(hmax);
289 fAmpLEDRec.AddAtAndExpand(grled,i);
083f0bec 290
291 delete [] xx;
292 delete [] yy;
293 delete [] xxled;
294 delete [] yyled;
295 delete prY;
296 delete prledY;
700e6b36 297
2e6a5ee0 298 // }
083f0bec 299 delete cfd;
300 delete qtcVScfd;
301 delete ledVScfd;
a2ad8166 302 }
700e6b36 303
a2ad8166 304}
305