Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DBCALClump.cc
Go to the documentation of this file.
1 /*
2  * DBCALClump.cc
3  *
4  * Created by Beni Zihlmann Tue Mar 12 2013
5  * version 0.1
6  *
7  */
8 
9 #include "BCAL/DBCALClump.h"
10 
11 #include <TF1.h>
12 #include <TGraph.h>
13 #include <TGraphErrors.h>
14 #include <TMath.h>
15 
16 DBCALClump::DBCALClump(vector<const DBCALHit*> U, vector<const DBCALHit*> D){
17 
18  HitsU = U;
19  HitsD = D;
20 
21 }
22 
24 
25  // now do some work on the Clumps up stream and down stream
26  // create clump profiles and determine Energy Sums
27 
28  //const DBCALHit *BcalMatrixU[48*4][4]; // Up stream
29  //const DBCALHit *BcalMatrixD[48*4][4]; // Down stream
30  float EU[48*4];
31  float ED[48*4];
32 
33  fillArrays(EU,ED);
34 
35  float EmaxU=0.;
36  float EmaxD=0.;
37 
38  // find Maximum of energy in Clump array
39  int MaxU=999;
40  int MaxD=999;
41  for (int i=0;i<48*4;i++){
42  if ((EU[i]>0) && (ED[i]>0)) {
43  if (EU[i]>EmaxU){
44  EmaxU = EU[i];
45  MaxU = i;
46  }
47  if (ED[i]>EmaxD){
48  EmaxD = ED[i];
49  MaxD = i;
50  }
51  }
52  }
53 
54  //cout<<MaxU<<" / "<<MaxD<<endl;
55  // use largest energy to define central bin
56  if (MaxU>MaxD){
57  MaxD = MaxU;
58  } else {
59  MaxU = MaxD;
60  }
61 
62  resetProfiles();
63 
64  int UCounter=0;
65  int DCounter=0;
66 
67  // set up a shower profile such that the maximum energy of the shower
68  // ends up in bin 20;
69  int CenterOffsetU = 20-MaxU;
70  int CenterOffsetD = 20-MaxD;
71 
72  // create shower profiles for up stream and down stream
73  int n=MaxU;
74  // start from center upwards
75  while (EU[n]>0.){
76  int idx = n+CenterOffsetU;
77  if (idx<0){
78  idx += MaxU;
79  }
80  //cout<<"U> "<<n<<" "<<EU[n]<<" "<<idx<<endl;
81  if ((idx>=60)||(idx<0)){
82  break;
83  }
84  ProfileU[idx] = EU[n]; // limit shower upper limit
85  UCounter++;
86  // section to get average mean time and average time difference for sector
87  int loc=999;
88  for (unsigned int k=0;k<Sector.size();k++){
89  if (Sector[k] == n){
90  loc = k;
91  break;
92  }
93  }
94  if (loc<999){
95  float Counter=0.;
96  while (Sector[loc]==n){
97  ProfileMT[idx] += MeanTime[loc];
98  ProfileTD[idx] += DeltaTime[loc];
99  loc++;
100  Counter += 1.0;
101  }
102  ProfileMT[idx] /= Counter;
103  ProfileTD[idx] /= Counter;
104  }
105 
106  n++;
107  if (n>=48*4){
108  n -= 48*4;
109  }
110  }
111  // now move from center down wards
112  n = MaxU-1;
113  if (n<0){
114  n += 48*4;
115  }
116  while (EU[n]>0.){
117  int idx = n+CenterOffsetU;
118  //cout<<"U< "<<n<<" "<<EU[n]<<" "<<idx<<endl;
119  if (idx>20){
120  idx -= (48*4);
121  }
122  if (idx<0){
123  break;
124  }
125  ProfileU[idx] = EU[n];
126  UCounter++;
127  // section to get average mean time and average time difference for sector
128  int loc=999;
129  for (unsigned int k=0;k<Sector.size();k++){
130  if (Sector[k] == n){
131  loc = k;
132  break;
133  }
134  }
135  if (loc<999){
136  float Counter=0.;
137  while (Sector[loc]==n){
138  ProfileMT[idx] += MeanTime[loc];
139  ProfileTD[idx] += DeltaTime[loc];
140  loc++;
141  Counter += 1.0;
142  }
143  ProfileMT[idx] /= Counter;
144  ProfileTD[idx] /= Counter;
145  }
146 
147  n--;
148  if (n<0){
149  n += 48*4;
150  }
151  }
152  n=MaxD;
153  // start from center upwards
154  while (ED[n]>0.){
155 
156  int idx = n+CenterOffsetD;
157  //cout<<"D> "<<n<<" "<<ED[n]<<" "<<idx<<endl;
158  if (idx<0){
159  idx += MaxD;
160  }
161  if ((idx>=60)||(idx<0)){
162  break;
163  }
164  ProfileD[idx] = ED[n];
165  DCounter++;
166  n++;
167  if (n>=48*4){
168  n -= 48*4;
169  }
170  }
171  // now move from center down wards
172  n = MaxD-1;
173  while (ED[n]>0.){
174 
175  int idx = n+CenterOffsetD;
176  //cout<<"D< "<<n<<" "<<ED[n]<<" "<<idx<<endl;
177  if (idx>20){
178  idx -= (48*4);
179  }
180  if (idx<0){
181  break;
182  }
183  ProfileD[idx] = ED[n];
184  DCounter++;
185  n--;
186  if (n<0){
187  n += 48*4;
188  }
189  }
190 
191  // at this point we have a shower profile for the upstream and down stream Clump
192 
193  if (0){
194 
195  for (int k=0; k<60; k++){
196  cout<<ProfileU[k]<<" ";
197  }
198  cout<<endl;
199  for (int k=0; k<60; k++){
200  cout<<ProfileD[k]<<" ";
201  }
202  cout<<endl;
203 
204  }
205 
206 
207  // AT THIS POINT CODE SHOULD GO THAT TESTS IF THERE IS MORE THAN ONE
208  // PARTICLE-SHOWER IN THE CLUMP!!!!!!
209  // THIS IS NOT YET DONE IN VERSION 0.1 each Clump is just one shower
210 
211  double x[60];
212  for (int k=0;k<60;k++){
213  x[k] = (double) k;
214  }
215  int low = 0;
216  while (ProfileU[low]==0.){
217  low++;
218  }
219  int high = low;
220  while (ProfileU[high]!=0){
221  high++;
222  }
223  low--;
224  high++;
225 
226  //double chi2U = 999.;
227  double posU = 20.;
228  if (UCounter>1){
229  TGraph* grU = new TGraph(60,x,ProfileU);
230  //grU->Fit("gaus","Q","r",(Double_t)low,(Double_t)high);
231 
232  grU->Fit("gaus","Q","r",15.,25.);
233  TF1 *f1U = grU->GetFunction("gaus");
234 
235  if (f1U) {
236  //chi2U = f1U->GetChisquare();
237  posU = f1U->GetParameter(1);
238  } else {
239  posU = 20.;
240  //chi2U = 999.;
241  }
242  if (TMath::Abs(posU-20.)>2.){ // this shower has a complicated shape
243  f1U->SetParameter(1,20.);
244  f1U->FixParameter(2,1.0);
245  //grU->Fit("gaus","Q","r",(Double_t)low,(Double_t)high);
246  grU->Fit("gaus","Q","r",15.,25.);
247  f1U = grU->GetFunction("gaus");
248  posU = f1U->GetParameter(1);
249  if (TMath::Abs(posU-20.)>2.){ // still no better result
250  posU = 20.; // fix it.
251  }
252  }
253  } else if (UCounter){
254  posU = 20.;
255  //chi2U = 999;
256  } else {
257  cout<<"Error no data to fit U: this should never happen!!!"<<endl;
258  posU = 0;
259  //chi2U=999999;
260  }
261 
262  low = 0;
263  while (ProfileU[low]==0.){
264  low++;
265  }
266  high = low;
267  while (ProfileU[high]!=0){
268  high++;
269  }
270  low--;
271  high++;
272 
273  //double chi2D = 999.;
274  double posD = 20.;
275  if (DCounter>1){
276  TGraph* grD = new TGraph(60,x,ProfileD);
277  //grD->Fit("gaus","Q","r",(Double_t)low,(Double_t)high);
278  grD->Fit("gaus","Q","r",15.,25.);
279 
280  TF1 *f1D = grD->GetFunction("gaus");
281  if (f1D){
282  //chi2D = f1D->GetChisquare();
283  posD = f1D->GetParameter(1);
284  } else {
285  posD = 20.;
286  //chi2D = 999.;
287  }
288  if (TMath::Abs(posD-20.)>2.){ // this shower has a complicated shape
289  f1D->SetParameter(1,20.);
290  f1D->FixParameter(2,1.0);
291  //grD->Fit("gaus","Q","r",(Double_t)low,(Double_t)high);
292  grD->Fit("gaus","Q","r",15.,25.);
293  f1D = grD->GetFunction("gaus");
294  posD = f1D->GetParameter(1);
295  if (TMath::Abs(posD-20.)>2.){ // still no better result
296  posD = 20.; // fix it.
297  }
298  }
299  } else if (DCounter){
300  posD = 20.;
301  //chi2D = 999.;
302  } else {
303  cout<<"Error no data to fit D: this should never happen!!!"<<endl;
304  posU = 0;
305  //chi2U=999999;
306  }
307 
308 
309  //cout<<chi2U<<" "<<chi2D<<" ";
310  //cout<<posU<<" "<<posD<<endl;
311 
312  //double UD_asymmetry = (posD-posU)/posD; // Note this quantity may be used for multi shower searches.
313 
314  // The following assumes that the first BCAL module (0) starts at 9 o'clock
315  // looking down stream and the counting is goind clock wise.
316  float phiD = (posD-20.+MaxD+0.5) * 2.*3.1415926/(48.*4.) ;
317  float phiU = (posU-20.+MaxU+0.5) * 2.*3.1415926/(48.*4.) ;
318  ClumpPhi.push_back((phiD+phiU)/2.);
319  //cout<<phid*180./3.1415926<<endl;
320 
321 
322  // calculate total energy of the shower from up and down stream clump energy
323  // using the geometri mean approach
324  double ESumU = 0.0;
325  double ESumD = 0.0;
326  for (unsigned int k=0;k<60;k++){
327  ESumU += ProfileU[k];
328  ESumD += ProfileD[k];
329  }
330 
331  // ATTENTION HARD CODED VALUES 195/300 = BcalLength/2/ATTENUATIONLENGTH !!!!!!!
332  double E = TMath::Sqrt(ESumU*ESumD)/TMath::Exp(-195./300.);
333  ClumpE.push_back(E);
334  // cout<<E<<endl;
335 
336  float MTaverage=0;
337  double Position[4] = {0.,0.,0.,0.};
338  double err[4] = {200.,200.,200.,200.};
339  double ex[4] = {0.,0.,0.,0.};
340  double Pcnt[4] = {0.,0.,0.,0.};
341  for (unsigned int k=0;k<MeanTime.size();k++){
342  MTaverage += MeanTime[k];
343  int idx = Layer[k];
344  Position[idx] += DeltaTime[k];
345  Pcnt[idx] += 1.;
346  }
347  if (MeanTime.size()){
348  MTaverage /= (float)MeanTime.size();
349  } else{
350  MTaverage = 0.;
351  }
352  ClumpMT.push_back(MTaverage);
353 
354  low=4;
355  high=0;
356  for (int k=0;k<4;k++){
357  if (Pcnt[k]>0.){
358  Position[k] /= Pcnt[k];
359  err[k] = 2.5; // Hard coded error
360  if (k<low){
361  low = k;
362  }
363  if (k>high){
364  high = k;
365  }
366  }
367  }
368 
369  double Pos=0;
370  double cnt = 0;
371  // don't do a fit with only two or less points
372  // fit approach seem to be too instable
373  // do at the moment do only average
374  if ((high-low)<6){
375  for (int k=0;k<4;k++){
376  if (Pcnt[k]>0){
377  Pos += Position[k];
378  cnt += 1.;
379  }
380  }
381  if (cnt>0){
382  Pos /= cnt;
383  }
384 
385  //cout<<" "<<Pos<<endl;
386  //for (int i=0;i<4;i++){
387  // cout<<Position[i]<<" ";
388  //}
389  //cout<<endl;
390 
391  } else { // do a fit with 3 or more points
392  double xx[4] = {1.,3.,8.,15.};
393  TGraphErrors* grT = new TGraphErrors(4,xx,Position,ex,err);
394  grT->Fit("pol1","Q","R",low,high);
395  TF1 * lfit = grT->GetFunction("pol1");
396  double offset = lfit->GetParameter(0);
397  //double slope = lfit->GetParameter(1);
398  //cout<<slope<<" * x + "<<offset<<endl;
399  //for (int i=0;i<4;i++){
400  // cout<<Position[i]<<" ";
401  // }
402  //cout<<endl;
403  Pos = offset;
404  }
405 
406  // adjust postion to be zero at the up stream end
407  // and full length at the down stream end
408  Pos = 195.-Pos; // hard code half length of module
409  ClumpPos.push_back(Pos);
410  //cout<<Pos<<endl;
411 
412 }
413 
414 void DBCALClump::fillArrays(float* EU, float* ED){
415 
416  // reset array
417  for (int i=0;i<48*4;i++){
418  EU[i] = 0.;
419  ED[i] = 0.;
420  }
421 
422  for (unsigned int i=0;i<HitsU.size();i++){
423  const DBCALHit *hit = HitsU[i];
424  int idx = (hit->sector-1) + (hit->module-1)*4;
425  EU[idx] += hit->E;
426  }
427  for (unsigned int i=0;i<HitsD.size();i++){
428  const DBCALHit *hit = HitsD[i];
429  int idx = (hit->sector-1) + (hit->module-1)*4;
430  ED[idx] += hit->E;
431  }
432 }
433 
435 
436  for (int k=0;k<60;k++){
437  ProfileU[k] = 0.0;
438  ProfileD[k] = 0.0;
439  ProfileMT[k] = 0.0;
440  ProfileTD[k] = 0.0;
441  }
442 
443 }
444 
vector< float > DeltaTime
Definition: DBCALClump.h:33
vector< float > ClumpPhi
Definition: DBCALClump.h:46
void fillArrays(float *, float *)
Definition: DBCALClump.cc:414
float E
Definition: DBCALHit.h:30
double ProfileMT[60]
Definition: DBCALClump.h:38
int module
Definition: DBCALHit.h:25
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
DBCALClump(vector< const DBCALHit * >, vector< const DBCALHit * >)
Definition: DBCALClump.cc:16
vector< float > ClumpPos
Definition: DBCALClump.h:45
vector< int > Sector
Definition: DBCALClump.h:34
Double_t ex[NCHANNELS]
vector< float > ClumpE
Definition: DBCALClump.h:43
double ProfileU[60]
Definition: DBCALClump.h:36
int sector
Definition: DBCALHit.h:27
void resetProfiles(void)
Definition: DBCALClump.cc:434
vector< int > Layer
Definition: DBCALClump.h:35
vector< const DBCALHit * > HitsU
Definition: DBCALClump.h:30
vector< float > ClumpMT
Definition: DBCALClump.h:44
vector< float > MeanTime
Definition: DBCALClump.h:32
double ProfileD[60]
Definition: DBCALClump.h:37
vector< const DBCALHit * > HitsD
Definition: DBCALClump.h:31
void AnalyzeClump()
Definition: DBCALClump.cc:23
double ProfileTD[60]
Definition: DBCALClump.h:39