AMSDST
myStatus.cxx
Go to the documentation of this file.
1 // Authors: M.Duranti, D.D'Urso - INFN di Perugia
2 #include "myStatus.h"
3 #include "myEvent.h"
4 #include "debug.h"
5 #include "TClass.h"
6 #include <iostream>
7 
8 using namespace std;
9 
10 //--------------------------------------------------------------------
11 
12 ClassImp(myStatus);
13 #ifdef _WITHGBATCH_
14 ClassImp(myStatusFiller);
15 #endif //#ifdef _WITHGBATCH_
16 
17 //--------------------------------------------------------------------
18 //#define PDEBUG
19 //#ifdef PDEBUG
20 //#define PRINTDEBUG printf("%s) This is the line number %d\n", __FILE__, __LINE__);
21 //#endif
22 
23 myStatus::myStatus(){
24 #ifdef PDEBUG
25  printf("In myStatus::myStatus\n");
26 #endif
27  PRINTDEBUG;
28  init();
29  PRINTDEBUG;
30 }
31 
32 myStatus::~myStatus(){
33 #ifdef PDEBUG
34  printf("In myStatus::~myStatus\n");
35 #endif
36  PRINTDEBUG;
37  Clear();
38 }
39 
40 void myStatus::Clear(Option_t* option){
41 #ifdef PDEBUG
42  printf("In myStatus::Clear\n");
43 #endif
44  PRINTDEBUG;
45  fStatus=0;
46  PRINTDEBUG;
47  GeneralStatus='\0';
48  Rigidity='\0';
49  EcalInfo='\0';
50  LogLTrd='\0';
51  TrkCharge='\0';
52  QTof='\0';
53  Detector='\0';
54  PRINTDEBUG;
55  EntryNumber=-9999;
56  TreeNumber=255;
57  PRINTDEBUG;
58  return;
59 }
60 
61 void myStatus::init(){
62 #ifdef PDEBUG
63  printf("In myStatus::init\n");
64 #endif
65  PRINTDEBUG;
66  Clear();
67  PRINTDEBUG;
68  return;
69 }
70 
71 UShort_t myStatus::IsValidTime(){
72 #ifdef PDEBUG
73  printf("In myStatus::IsValidTime\n");
74 #endif
75 PRINTDEBUG;
76 
77  myRTI *myrti = myEvent::gethead()->RTI;
78  if(!myrti)return 0; //Time no information
79  if(myrti->good==0){
80 #ifdef PDEBUG
81  cout << " Abandon Second "<< endl;
82 #endif
83  return 0;//Abandon Second
84  }else
85  return 1;
86 }
87 
88 UShort_t myStatus::CodeRigidity(Double_t R, bool IsRigidity){
89 #ifdef PDEBUG
90  printf("In myStatus::CodeRigidity\n");
91 #endif
92  PRINTDEBUG;
93 
94  if(R==0) return 0;
95  int sign=(R>0)?1:0;
96  int return_value=-9999;
97  int power_limit=(IsRigidity==true)?7:8;
98  if(TMath::Abs(R) > pow(2.,power_limit)){
99  return_value=power_limit;
100  } else {
101  double log2r=(TMath::Abs(R)>=1)?log2(TMath::Abs(R)):0;
102  return_value=(int)log2r;
103  }
104  if(IsRigidity==true){
105  return_value=(return_value<<1);
106  return_value+=sign;
107  }
108 #ifdef PDEBUG
109  cout<<" return_value = "<<return_value<<endl;
110 #endif
111  return (UShort_t) return_value;
112 }
113 
114 UShort_t myStatus::CodeCharge(Double_t Q){
115 #ifdef PDEBUG
116  printf("In myStatus::CodeCharge\n");
117 #endif
118  PRINTDEBUG;
119  UShort_t return_value=0;
120  PRINTDEBUG;
121  double remapped= Q*3.75;
122  if(remapped>=15) remapped=15;
123  return_value= (UShort_t)remapped;
124  return return_value;
125 }
126 
127 UShort_t myStatus::CodeBDT(Double_t BDT){
128 #ifdef PDEBUG
129  printf("In myStatus::CodeBDT\n");
130 #endif
131  PRINTDEBUG;
132  UShort_t return_value=0;
133  PRINTDEBUG;
134  int sign=(BDT>0)?1:0;
135  double bdt_remapped= log2((1. - TMath::Abs(BDT))*128.);
136  return_value= (UShort_t)bdt_remapped;
137  return_value+=sign;
138  return return_value;
139 }
140 
141 UShort_t myStatus::CodeLikeL(Double_t LikeL){
142 #ifdef PDEBUG
143  printf("In myStatus::CodeLikeL\n");
144 #endif
145  PRINTDEBUG;
146  UShort_t return_value=0;
147  PRINTDEBUG;
148  if(LikeL<0) return return_value;
149  double remapped= LikeL*7.5;
150  if(remapped>=15) remapped=15;
151  return_value= (UShort_t)remapped;
152  return return_value;
153 }
154 
155 //-----------------------------------------------------------------------
156 
157 #ifdef _WITHGBATCH_
158 
159 myStatusFiller::myStatusFiller(){
160 #ifdef PDEBUG
161  printf("In myStatusFiller::myStatusFiller\n");
162 #endif
163  PRINTDEBUG;
164  init();
165  PRINTDEBUG;
166 }
167 
168 myStatusFiller::~myStatusFiller(){
169 #ifdef PDEBUG
170  printf("In myStatusFiller::~myStatusFiller\n");
171 #endif
172  PRINTDEBUG;
173  Clear();
174 }
175 
176 void myStatusFiller::Clear(Option_t* option){
177 #ifdef PDEBUG
178  printf("In myStatusFiller::Clear\n");
179 #endif
180  PRINTDEBUG;
181  return;
182 }
183 
184 void myStatusFiller::init(){
185 #ifdef PDEBUG
186  printf("In myStatusFiller::init\n");
187 #endif
188  PRINTDEBUG;
189  Clear();
190  PRINTDEBUG;
191  return;
192 }
193 
194 void myStatusFiller::Fill(int i_part_max_r, int i_part_max_qbh, int i_part_max_e, int i_part_max_r_trd){
195 #ifdef PDEBUG
196  printf("In myStatusFiller::Fill\n");
197 #endif
198  myEvent *mev=myEvent::gethead();
199  PRINTDEBUG;
200 #ifdef PDEBUG
201  cout<<" Event "<<mev->Event<<" Run "<<mev->Run<<endl;
202 #endif
203  int nParticle = mev->nParticle();
204  PRINTDEBUG;
205 
206  GeneralStatus=(UChar_t)IsValidTime();
207 
208  PRINTDEBUG;
209 
210  if(nParticle==0)
211  GeneralStatus=(0<<1);
212  else if(nParticle>=1 && nParticle<=2)
213  GeneralStatus=(nParticle<<1);
214  else
215  GeneralStatus=(3<<1);
216 
217 
218  PRINTDEBUG;
219  int Refit=1;
220  Double_t Rigidity_inn=0.;
221  Double_t Rigidity=0.;
222  Double_t InnerQ=0.;
223  Double_t QL1=0.;
224  if(i_part_max_r>=0){
225  myTrTrack *mtr=mev->pmyParticle(i_part_max_r)->pmyTrTrackByMySelf();
226  PRINTDEBUG;
227  if(mtr){
228  AMSEventR *ev=AMSEventR::Head();
229  TrTrackR *trk=ev->pTrTrack(mtr->i_myself);
230  if (trk) {
231  int itrtrack_inn=trk->iTrTrackPar(1, 3, Refit);
232  if(itrtrack_inn>=0){
233  Rigidity_inn=trk->GetRigidity(itrtrack_inn);
234  }
235 
236  int itrtrack=trk->iTrTrackPar(1, 0, Refit);
237  if(itrtrack>=0){
238  Rigidity=trk->GetRigidity(itrtrack);
239  }
240 
241 #ifdef PDEBUG
242  bool HasL1 = ((mtr->BitPatternJ)&(1<<0))>>0;
243  bool HasL2 = ((mtr->BitPatternJ)&(1<<1))>>1;
244  bool HasL9 = ((mtr->BitPatternJ)&(1<<8))>>8;
245  // if(HasL1&&HasL2&&HasL9 ) exit(1);
246 #endif
247 
248  InnerQ=mtr->InnerQ;
249  QL1=mtr->LayerJQ[0];
250  }
251  }
252  PRINTDEBUG;
253  StoreRigidity(Rigidity_inn,
254  Rigidity);
255 
256  TrkCharge=(UChar_t)CodeCharge(InnerQ);
257  UChar_t q_remapped=(UChar_t) CodeCharge(QL1);
258  TrkCharge+=(UChar_t) (q_remapped<<4);
259 #ifdef PDEBUG
260  cout<<" InnerQ = "<<InnerQ
261  <<" QL1 = "<<QL1
262  <<(int)TrkCharge<<endl;
263 #endif
264 
265  StoreDetectorAcceptance(i_part_max_r);
266 
267  }
268  PRINTDEBUG;
269 
270  if(i_part_max_qbh>=0){
271  myBetaHPlus *mybetahplus=mev->pmyParticle(i_part_max_qbh)->pmyBetaHPlusByMySelf();
272  PRINTDEBUG;
273  if(mybetahplus){
274  // AMSEventR *ev=AMSEventR::Head();
275  // BetaHR *betah=ev->pBetaH(mybetahplus->i_myself);
276  PRINTDEBUG;
277  // if(betah){
278  // QTof=(UChar_t) CodeCharge( mybetahplus->GetTofCharge(1));
279  // UChar_t remapped=(UChar_t) CodeCharge(mybetahplus->GetTofCharge(2));
280  // QTof += (remapped<<4);
281 #ifdef PDEBUG
282  // cout<<" QUTof = "<<mybetahplus->GetTofCharge(1)
283  // <<" QLTof = "<<mybetahplus->GetTofCharge(2)
284  // <<(int)QTof<<endl;
285 #endif
286  // }
287  PRINTDEBUG;
288  }
289  PRINTDEBUG;
290  }
291  PRINTDEBUG;
292 
293  if(i_part_max_e>=0 ){
294  myEcalShower *myecal=mev->pmyParticle(i_part_max_e)->pmyEcalShowerByMySelf();
295  myEcalShowerPlus *myecalplus=mev->pmyParticle(i_part_max_e)->pmyEcalShowerPlusByMySelf();
296  if (myecal && myecalplus) {
297  EcalInfo=(UChar_t) CodeRigidity(myecal->EnergyD, false);
298  UChar_t bdt_remapped=(UChar_t) CodeBDT(myecalplus->BDT);
299  EcalInfo += (bdt_remapped<<4);
300 #ifdef PDEBUG
301  cout<<" EnergyD = "<<myecal->EnergyD<<" "
302  <<" BDT = "<<myecalplus->BDT<<" "<<(int)EcalInfo<<endl;
303 #endif
304  }
305 
306  }
307 
308  PRINTDEBUG;
309 
310  if(i_part_max_r_trd>=0){
311  myTrdK *mytrdplus=mev->pmyParticle(i_part_max_r_trd)->pmyTrdKFromTrTrackByMySelf();
312  if (mytrdplus) {
313  LogLTrd=(UChar_t) CodeLikeL(mytrdplus->LHR[0]);
314  UChar_t remapped=(UChar_t) CodeLikeL(mytrdplus->LHR[1]);
315  LogLTrd += (remapped<<4);
316 #ifdef PDEBUG
317  cout<<" TrdKLHR[0] = "<<mytrdplus->LHR[0]<<" "
318  <<" TrdKLHR[1] = "<<mytrdplus->LHR[1]
319  <<(int)LogLTrd<<endl;
320 #endif
321  }
322  }
323 
324  return;
325 }
326 
327 void myStatusFiller::StoreRigidity(Double_t R_inn, Double_t R_max){
328 #ifdef PDEBUG
329  printf("In myStatusFiller::StoreRigidity\n");
330 #endif
331 PRINTDEBUG;
332 
333  UShort_t r_inn= CodeRigidity(R_inn);
334  UShort_t r_maxspan= CodeRigidity(R_max);
335  r_maxspan=(r_maxspan<<4);
336  Rigidity = r_inn;
337  Rigidity += r_maxspan;
338 #ifdef PDEBUG
339  cout<<" Rigidity = "<<r_inn<<" + ( "<<r_maxspan<<" ) = "<<(int) Rigidity<<endl;
340 #endif
341 }
342 
343 void myStatusFiller::StoreDetectorAcceptance(int i_part_max_r){
344 
345  //Bit Map: 0 L1, 1 TRD, 2 UTof, 3 L2, 4 LTof, 5 L9, 6 Ecal
346  float margin=0.;
347  if(i_part_max_r<0) return;
348  myEvent *mev=myEvent::gethead();
349  myTrTrack *mytrk=mev->pmyParticle(i_part_max_r)->pmyTrTrackByMySelf();
350  PRINTDEBUG;
351  if(mytrk){
352  AMSEventR *ev=AMSEventR::Head ();
353  TrTrackR *trk=ev->pTrTrack(mytrk->i_myself);
354  PRINTDEBUG;
355  if(trk){
356  AMSPoint point;
357  AMSDir dir;
358  int i_layer=1;
359 
360  trk->InterpolateLayerJ(i_layer,point,dir,0);
361  if (IsInsideLayer(i_layer,point,margin))
362  Detector=1;
363 
364  if(InTrdAcceptance(trk, mytrk->kDef))
365  Detector |= (1<<1);
366 
367  trk->Interpolate(66, point, dir, mytrk->kDef);
368  if( TMath::Abs(point.x())<65 && TMath::Abs(point.y())<65)
369  Detector |= (1<<2);
370 
371  i_layer=2;
372  trk->InterpolateLayerJ(i_layer,point,dir,0);
373  if (IsInsideLayer(i_layer,point,margin))
374  Detector |= (1<<3);
375 
376  trk->Interpolate(-66, point, dir, mytrk->kDef);
377  if( TMath::Abs(point.x())<65 && TMath::Abs(point.y())<65)
378  Detector |= (1<<4);
379 
380  i_layer=9;
381  trk->InterpolateLayerJ(9,point,dir,0);
382  if (IsInsideLayer(9,point,margin))
383  Detector |= (1<<5);
384 
385  AMSPoint APointTrackEntryEcal;
386  AMSDir ADirTrackEntryEcal;
387  trk->Interpolate(-142.79, APointTrackEntryEcal, ADirTrackEntryEcal, mytrk->kDef);
388 
389  AMSPoint APointTrackExitEcal;
390  AMSDir ADirTrackExitEcal;
391  trk->Interpolate(-160.00, APointTrackExitEcal, ADirTrackExitEcal, mytrk->kDef);
392 
393  if ( ( TMath::Abs(APointTrackEntryEcal.x())<35 && TMath::Abs(APointTrackEntryEcal.y())<35)
394  && ( TMath::Abs(APointTrackExitEcal.x())<35 && TMath::Abs(APointTrackExitEcal.y())<35)
395  )
396  Detector |= (1<<6);
397 
398  }
399 
400  }
401 }
402 
403 
404 bool myStatusFiller::IsInsideLayer(int jlayer, AMSPoint point, float margin) {
405 
406  // float tracker_layers_z[9] = {158.920,53.060,29.228,25.212,1.698,-2.318,-25.212,-29.228,-135.882};
407  float tracker_planes_edges[9][4] = {
408  {-62.14, -47.40, 62.14, 47.40},
409  {-62.14, -40.10, 62.14, 40.10},
410  {-49.70, -43.75, 49.70, 43.75},
411  {-49.72, -43.75, 49.72, 43.75},
412  {-49.71, -36.45, 49.70, 36.45},
413  {-49.72, -36.45, 49.72, 36.45},
414  {-49.72, -43.75, 49.71, 43.75},
415  {-49.72, -43.75, 49.71, 43.75},
416  {-45.62, -29.48, 45.55, 29.53}
417  };
418 
419  bool isinlayer = false;
420  int ilayer = jlayer - 1;
421  if ( (point.x()>tracker_planes_edges[ilayer][0]+margin)&&
422  (point.x()<tracker_planes_edges[ilayer][2]-margin)&&
423  (point.y()>tracker_planes_edges[ilayer][1]+margin)&&
424  (point.y()<tracker_planes_edges[ilayer][3]-margin) ) {
425  if ((ilayer+1)==9) isinlayer = true;
426  else {
427  // circle
428  if ( (sqrt(point.x()*point.x()+point.y()*point.y())<
429  tracker_planes_edges[ilayer][2]-margin) )
430  isinlayer = true;
431  }
432  }
433  return isinlayer;
434 }
435 
436 
437 bool myStatusFiller::InTrdAcceptance(TrTrackR* this_trk, int id_maxspan){
438  //check if particle passing inside acceptance of the 20 layers TRD
439 
440  if(!this_trk) return true;//return true to have the cutnotpassed evaluated more reasonably
441  //only top layer and bottom layer checked
442  //this is a gross check. A more fine tuned check is going to be provided by Aachen group
443 
444  int nTrdBottom = 9;
445  float AccepBottomX[] = {+40, +78, +78, +40, -40, -78, -78, -40, +40};
446  float AccepBottomY[] = {+76, +35, -35, -76, -76, -35, +35, +76, +76};
447  float TrdBottomZ = 86.725;
448 
449  // int nTrdCenter = 9;
450  // float AccepCenterX[] = {-80.0, -47.0, 47.0, 80.0, 80.0, 47.0, -47.0, -80.0, -80.0};
451  // float AccepCenterY[] = { 43.5, 75.5, 75.5, 43.5, -43.5, -75.5, -75.5, -43.5, 43.5};
452  float TrdCenterZ = 0.5 * (141.825 + 86.725);
453 
454  int nTrdTop = 37;
455  float AccepTopX[] = {-99.0,-89.0,-89.0,-78.7,-78.7,-67.8,-67.8,-57.7,-57.7, 57.7, 57.7, 67.8, 67.8, 78.7, 78.7, 89.0, 89.0, 99.0,
456  99.0, 89.0, 89.0, 78.7, 78.7, 67.8, 67.8, 57.7, 57.7,-57.7,-57.7,-67.8,-67.8,-78.7,-78.7,-89.0,-89.0,-99.0,-99.0};
457  float AccepTopY[] = { 54.5, 54.5, 62.5, 62.5, 74.0, 74.0, 84.0, 84.0, 95.3, 95.3, 84.0, 84.0, 74.0, 74.0, 62.5, 62.5, 54.5, 54.5,
458  -51.7,-51.7,-62.2,-62.2,-72.0,-72.0,-82.5,-82.5,-92.5,-92.5,-82.5,-82.5,-72.0,-72.0,-62.2,-62.2,-51.7,-51.7, 54.5};
459  float TrdTopZ = 141.825;
460 
461  AMSPoint pTop, pCenter, pBottom;
462  AMSDir dTop, dCenter, dBottom;
463 
464  this_trk->Interpolate(TrdBottomZ, pBottom, dBottom, id_maxspan);
465  this_trk->Interpolate(TrdCenterZ, pCenter, dCenter, id_maxspan);
466  this_trk->Interpolate(TrdTopZ, pTop, dTop, id_maxspan);
467 
468  bool passTrdBottom = TMath::IsInside((float)pBottom.x(), (float)pBottom.y(), nTrdBottom, AccepBottomX, AccepBottomY);
469  // bool passTrdCenter = TMath::IsInside((float)pCenter.x(), (float)pCenter.y(), nTrdCenter, AccepCenterX, AccepCenterY);
470  bool passTrdTop = TMath::IsInside((float)pTop.x(), (float)pTop.y(), nTrdTop, AccepTopX, AccepTopY);
471 
472  return (passTrdTop && passTrdBottom);
473 }
474 
475 #endif //#ifdef _WITHGBATCH_
Generated on Tue Dec 29 2015 23:02:50 for AMSDST by   doxygen 1.8.3