Line data Source code
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 : //
17 : // TRD dEdx recon utils
18 : // xx
19 : // xx
20 : // xx
21 : // xx
22 : //
23 : // Xianguo Lu
24 : // lu@physi.uni-heidelberg.de
25 : // Xianguo.Lu@cern.ch
26 : //
27 : //
28 :
29 : #include "TF1.h"
30 : #include "TFile.h"
31 : #include "TH1D.h"
32 : #include "TH2D.h"
33 : #include "THnSparse.h"
34 : #include "TMath.h"
35 : #include "TMatrixD.h"
36 : #include "TMinuit.h"
37 : #include "TObjArray.h"
38 : #include "TRandom3.h"
39 : #include "TStopwatch.h"
40 : #include "TVectorD.h"
41 :
42 : #include "TTreeStream.h"
43 :
44 : #include "AliCDBId.h"
45 : #include "AliCDBMetaData.h"
46 : #include "AliCDBStorage.h"
47 : #include "AliESDEvent.h"
48 : #include "AliESDfriendTrack.h"
49 : #include "AliESDtrack.h"
50 : #include "AliTRDcalibDB.h"
51 : #include "AliTRDCalROC.h"
52 : #include "AliTRDtrackV1.h"
53 :
54 : #include "AliTRDdEdxBaseUtils.h"
55 : #include "AliTRDdEdxReconUtils.h"
56 :
57 : #define EPSILON 1e-8
58 :
59 : Int_t AliTRDdEdxReconUtils::ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
60 : {
61 : //
62 : //apply calibration on arrayQ
63 : //
64 208 : if(!cobj){ printf("AliTRDdEdxReconUtils::ApplyCalib error gain array null!!\n"); exit(1);}
65 :
66 104 : TVectorD tmpq(arrayQ->GetNrows());
67 104 : TVectorD tmpx(arrayX->GetNrows());
68 : Int_t ncls = 0;
69 :
70 208 : const TVectorD * gain = (TVectorD*) cobj->At(0);
71 9096 : for(Int_t ii=0; ii<nc0; ii++){
72 8888 : const Double_t dq = (*arrayQ)[ii];
73 8888 : const Int_t xx = (Int_t)(*arrayX)[ii];
74 8888 : const Double_t gg = (*gain)[xx];
75 :
76 4444 : if(gg<EPSILON){
77 0 : continue;
78 : }
79 :
80 8888 : tmpq[ncls] = dq*gg;
81 8888 : tmpx[ncls] = xx;
82 4444 : ncls++;
83 4444 : }
84 :
85 104 : (*arrayQ)=tmpq;
86 104 : (*arrayX)=tmpx;
87 :
88 : return ncls;
89 104 : }
90 :
91 : Double_t AliTRDdEdxReconUtils::ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
92 : {
93 : //
94 : //template for cookdedx
95 : //
96 104 : if(cobj){
97 104 : if(arrayQ && arrayX){
98 104 : ncluster = ApplyCalib(ncluster, arrayQ, arrayX, cobj);
99 : }
100 : else{
101 0 : printf("AliTRDdEdxReconUtils::ToyCook arrayQ arrayX null, applycalib can not be applied!\n"); exit(1);
102 : }
103 104 : }
104 :
105 : Double_t lowFrac =-999, highFrac = -999;
106 104 : if(kinvq){
107 0 : lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
108 0 : }
109 : else{
110 104 : lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
111 : }
112 :
113 104 : Double_t meanQ = AliTRDdEdxBaseUtils::TruncatedMean(ncluster, arrayQ->GetMatrixArray(), lowFrac, highFrac);
114 104 : if(kinvq){
115 104 : if(meanQ>EPSILON){
116 0 : meanQ = 1/meanQ;
117 0 : }
118 : }
119 :
120 104 : return meanQ;
121 : }
122 :
123 : Double_t AliTRDdEdxReconUtils::CombineddEdx(const Bool_t kinvq, Int_t &concls, TVectorD *coarrayQ, TVectorD *coarrayX, const Int_t tpcncls, const TVectorD *tpcarrayQ, const TVectorD *tpcarrayX, const Int_t trdncls, const TVectorD *trdarrayQ, const TVectorD *trdarrayX)
124 : {
125 : //
126 : //combine tpc and trd dedx
127 : //
128 :
129 0 : for(Int_t iq=0; iq<tpcncls; iq++){
130 0 : (*coarrayQ)[iq]=(*tpcarrayQ)[iq];
131 0 : if(tpcarrayX && trdarrayX && coarrayX){
132 0 : (*coarrayX)[iq]=(*tpcarrayX)[iq];
133 0 : }
134 : }
135 0 : for(Int_t iq=0; iq<trdncls; iq++){
136 0 : (*coarrayQ)[tpcncls+iq]=(*trdarrayQ)[iq];
137 0 : if(tpcarrayX && trdarrayX && coarrayX){
138 0 : (*coarrayX)[tpcncls+iq]=159+(*trdarrayX)[iq];
139 0 : }
140 : }
141 :
142 0 : concls=trdncls+tpcncls;
143 :
144 0 : const Double_t coQ = ToyCook(kinvq, concls, coarrayQ, coarrayX);
145 :
146 0 : return coQ;
147 : }
148 :
149 : Double_t AliTRDdEdxReconUtils::GetPadGain(const Int_t det, const Int_t icol, const Int_t irow)
150 : {
151 : //
152 : //get pad calibration
153 : //
154 99864 : AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
155 49932 : if (!calibration) {
156 0 : printf("AliTRDdEdxReconUtils::GetPadCalib No AliTRDcalibDB instance available\n"); exit(1);
157 : }
158 49932 : AliTRDCalROC * calGainFactorROC = calibration->GetGainFactorROC(det);
159 49932 : if(!calGainFactorROC){
160 0 : printf("AliTRDdEdxReconUtils::GetPadCalib no calGainFactorROC!\n"); exit(1);
161 : }
162 :
163 : Double_t padgain = -999;
164 99864 : if( icol >= 0 &&
165 99864 : icol < calGainFactorROC->GetNcols() &&
166 49932 : irow >=0 &&
167 49932 : irow < calGainFactorROC->GetNrows()){
168 49932 : padgain = calGainFactorROC->GetValue(icol, irow);
169 49932 : if(padgain<EPSILON){
170 0 : printf("AliTRDdEdxReconUtils::GetPadGain padgain 0! %f %f -- %d %d %d -- %d %d\n", padgain, EPSILON, det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows()); exit(1);
171 : }
172 : }
173 : else{
174 : //printf("\nAliTRDdEdxReconUtils::GetPadGain warning!! indices out of range %d %d %d -- %d %d\n\n", det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows() );
175 : }
176 :
177 49932 : return padgain;
178 : }
179 :
180 : Double_t AliTRDdEdxReconUtils::GetRNDClusterQ(AliTRDcluster *cl, const Double_t baseline)
181 : {
182 : //
183 : //get cluter q from GetRawQ, apply baseline and Kr pad-calibration
184 : //
185 :
186 17920 : const Int_t det = cl->GetDetector();
187 8960 : const Int_t pad3col = cl->GetPadCol();
188 8960 : const Int_t padrow = cl->GetPadRow();
189 :
190 : Double_t rndqsum = 0;
191 143360 : for(Int_t ii=0; ii<7; ii++){
192 62720 : if(cl->GetSignals()[ii] < EPSILON){//bad pad marked by electronics
193 : continue;
194 : }
195 :
196 49932 : const Int_t icol = pad3col+(ii-3);
197 49932 : const Double_t padgain = GetPadGain(det, icol, padrow);
198 49932 : if(padgain<0){//indices out of range, pad3col near boundary case
199 0 : continue;
200 : }
201 :
202 49932 : const Double_t rndsignal = (cl->GetSignals()[ii] - baseline )/(AliTRDdEdxBaseUtils::IsPadGainOn()? padgain : 1);
203 :
204 : //sum it anyway even if signal below baseline, as long as the total is positive
205 49932 : rndqsum += rndsignal;
206 49932 : }
207 :
208 8960 : return rndqsum;
209 : }
210 :
211 : Double_t AliTRDdEdxReconUtils::GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb)
212 : {
213 : //
214 : //get cluster charge
215 : //
216 : Double_t dq = 0;
217 : AliTRDcluster *cl = 0x0;
218 :
219 : const Double_t baseline = 10;
220 :
221 : //GetRNDClusterQ(cl)>0 ensures that the total sum of q is above baseline*NsignalPhysical.
222 15018 : cl = seed->GetClusters(itb); if(cl && GetRNDClusterQ(cl, baseline)>0 ) dq+= GetRNDClusterQ(cl, baseline);
223 6714 : cl = seed->GetClusters(itb+AliTRDseedV1::kNtb); if(cl && GetRNDClusterQ(cl, baseline)>0 ) dq+= GetRNDClusterQ(cl, baseline);
224 :
225 6386 : dq /= AliTRDdEdxBaseUtils::Getdldx(seed);
226 :
227 6386 : dq /= AliTRDdEdxBaseUtils::QScale();
228 :
229 12772 : if(kinvq){
230 6386 : if(dq>EPSILON){
231 0 : dq = 1/dq;
232 0 : }
233 : }
234 :
235 6386 : return dq;
236 : }
237 :
238 : Int_t AliTRDdEdxReconUtils::GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0, Int_t timeBin1, Int_t tstep)
239 : {
240 : //
241 : //return nclustter
242 : //(if kinvq, return 1/q array), size of array must be larger than 31*6
243 : //
244 208 : if(!arrayQ || arrayQ->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
245 0 : printf("AliTRDdEdxReconUtils::GetArrayClusterQ error arrayQ null or size too small! %d\n", arrayQ? arrayQ->GetNrows() : -999); exit(1);
246 : }
247 208 : if(!arrayX || arrayX->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
248 0 : printf("AliTRDdEdxReconUtils::GetArrayClusterQ error arrayX null or size too small! %d\n", arrayX? arrayX->GetNrows() : -999); exit(1);
249 : }
250 :
251 : const Int_t mintb = 0;
252 : const Int_t maxtb = AliTRDseedV1::kNtb-1;
253 208 : if(timeBin0<mintb) timeBin0=mintb;
254 208 : if(timeBin1>maxtb) timeBin1=maxtb;
255 104 : if(tstep<=0) tstep=1;
256 :
257 : //============
258 : Int_t tbN=0;
259 104 : Double_t tbQ[200];
260 104 : Int_t tbBin[200];
261 :
262 1456 : for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
263 624 : const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
264 624 : if(!seed)
265 418 : continue;
266 :
267 206 : const Int_t det = seed->GetDetector();
268 :
269 13184 : for(Int_t itb=timeBin0; itb<=timeBin1; itb+=tstep){
270 6386 : const Double_t dq = GetClusterQ(kinvq, seed, itb);
271 6386 : if(dq<EPSILON)
272 1942 : continue;
273 :
274 4444 : const Int_t gtb = det * AliTRDseedV1::kNtb + itb;
275 :
276 4444 : tbQ[tbN]=dq;
277 4444 : tbBin[tbN]=gtb;
278 4444 : tbN++;
279 4444 : }
280 206 : }
281 :
282 : Int_t ncls = 0;
283 9096 : for(Int_t iq=0; iq<tbN; iq++){
284 4444 : if(tbQ[iq]<EPSILON)
285 : continue;
286 :
287 4444 : (*arrayQ)[ncls] = tbQ[iq];
288 4444 : (*arrayX)[ncls] = tbBin[iq];
289 :
290 4444 : ncls++;
291 4444 : }
292 :
293 : static Int_t kprint = 100;
294 104 : if(kprint<0){
295 0 : printf("\nAliTRDdEdxReconUtils::GetArrayClusterQ raw cluster-Q\n");
296 0 : for(Int_t iq=0; iq<ncls; iq++){
297 0 : const Int_t ichamber = AliTRDdEdxBaseUtils::ToLayer((*arrayX)[iq]);
298 0 : const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
299 0 : if(!seed){
300 0 : printf("error seed null!!\n"); exit(1);
301 : }
302 0 : const Double_t rawq = (*arrayQ)[iq] * 45. * AliTRDdEdxBaseUtils::Getdldx(seed);
303 0 : printf("esdid=%d; chamber=%d; timebin=%d; rawq= %.3f; myq[%d]= %e;\n", trdtrack->GetESDid(), ichamber, AliTRDdEdxBaseUtils::ToTimeBin((*arrayX)[iq]), rawq, iq, (*arrayQ)[iq]);
304 : }
305 0 : printf("\n");
306 0 : }
307 104 : kprint++;
308 :
309 104 : return ncls;
310 104 : }
311 :
312 : Int_t AliTRDdEdxReconUtils::UpdateArrayX(const Int_t ncls, TVectorD* arrayX)
313 : {
314 : //
315 : //arrayX det*Ntb+itb -> itb
316 : //
317 :
318 208 : TVectorD countChamber(6);
319 9096 : for(Int_t ii = 0; ii<ncls; ii++){
320 8888 : const Int_t xx = (Int_t)(*arrayX)[ii];
321 4444 : const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(xx);
322 :
323 8888 : const Double_t ich = AliTRDgeometry::GetLayer(idet);
324 8888 : const Double_t itb = AliTRDdEdxBaseUtils::ToTimeBin(xx);
325 8888 : (*arrayX)[ii] = ich*AliTRDseedV1::kNtb+itb;
326 :
327 8888 : countChamber[ich] = 1;
328 : }
329 :
330 104 : const Double_t nch = countChamber.Sum();
331 104 : return (Int_t) nch;
332 104 : }
333 :
334 :
|
