AMSDST
myRichRingPlus.cxx
Go to the documentation of this file.
1 // Authors: M.Duranti - INFN di Perugia
2 #include "myRichRingPlus.h"
3 #include "debug.h"
4 #include "TClass.h"
5 
6 using namespace std;
7 
8 //--------------------------------------------------------------------
9 
10 ClassImp(myRichRingPlus);
11 
12 //---------------------------------------------------------------------
13 
14 // FEW VARIABLES FOR GEOMTRICAL SELECTION ON RICH ------
15 const int grid_side_length=11;
16 const int n_tiles=grid_side_length*grid_side_length;
17 const Double_t tile_width=0.1+11.5;
18 // -----------------------------------------------------
19 
20 myRichRingPlus::myRichRingPlus(){
21 #ifdef PDEBUG
22  printf("In myRichRingPlus::myRichRingPlus\n");
23 #endif
24  PRINTDEBUG;
25  init();
26  PRINTDEBUG;
27 }
28 
29 myRichRingPlus::~myRichRingPlus(){
30 #ifdef PDEBUG
31  printf("In myRichRingPlus::~myRichRingPlus\n");
32 #endif
33  PRINTDEBUG;
34 }
35 
36 void myRichRingPlus::Clear(Option_t* option){
37 #ifdef PDEBUG
38  printf("In myRichRingPlus::Clear\n");
39 #endif
40  PRINTDEBUG;
41  Status = 999999;
42  Beta = -999999;
43  Theta = -999999;
44  JustOneRichParticle = false;
45  TheCorrectTrack = false;
46  _IsClean = false;
47  FiducialVolumeLoose = false;
48  FiducialVolumeTight = false;
49  EnoughUsedHitsInRing = false;
50  NoMirroredPhotons = false;
51  EnoughExpetedHitsInRing = false;
52  GoodRingProbability = false;
53  RingWidthSmallEnough = false;
54  EnoghPMTsInEvent = false;
55  EnoughPMTsInEventForRing = false;
56  EnoughPhotoElectronsWRTCollected = false;
57  BetaDiscrepancy = false;
58 
59  PRINTDEBUG;
60 
61  Used = -999999;
62  UsedM = -999999;
63  Prob = -999999;
64  fill_n(AMSTrPars, 5, -999999);
65  fill_n(TrRadPos, 3, -999999);
66  fill_n(TrPMTPos, 3, -999999);
67  BetaConsistency = -999999;
68  _betaCorrection = -999999;
69  PhotoElectrons = -999999;
70  ExpectedPhotoElectrons = -999999;
71  PMTs = -999999;
72  PMTChargeConsistency = -999999;
73  _PmtCorrectionsFailed = -999999;
74  TileIndex = -999999;
75  CollectedPhotoElectrons = -999999;
76  UsedHits = -999999;
77  ChargeCorrections=false;
78  _RingWidth = -999999;
79  NpExp = -999999;
80  NpCol = -999999;
81  BetaExpectedResolution = -999999;
82  BetaExpectedRms = -999999;
83  ChargeExpectedResolution = -999999;
84  ChargeExpectedRms = -999999;
85 
86  kind = -999999;
87  fill_n(mindist, 2, -999999);
88  fill_n(empoint, 3, -999999);
89  fill_n(emdir, 2, -999999);
90  fill_n(nphoel,2, -999999);
91  fill_n(betam,2, -999999);
92  fill_n(npmthits,2, -999999);
93  fill_n(npmts,2, -999999);
94 
95  PRINTDEBUG;
96 
97  return;
98 }
99 
100 void myRichRingPlus::init(){
101 #ifdef PDEBUG
102  printf("In myRichRingPlus::init\n");
103 #endif
104  PRINTDEBUG;
105  Clear();
106  PRINTDEBUG;
107  return;
108 }
109 
110 //---------------------------------------------------------------------
111 
112 #ifdef _WITHGBATCH_
113 
114 //---------------------------------------------------------------------
115 
116 Int_t richrec(TrTrackR const*, AMSPoint&, AMSDir& , AMSPoint & , double* , int* , int* , double* );
117 
118 //---------------------------------------------------------------------
119 
120 ClassImp(myRichRingPlusFiller);
121 
122 //--------------------------------------------------------------------
123 
124 static Double_t betaWidth(RichRingR* ring);
125 static int getPMTs(bool countCrossed=true);
126 static int getPMTsForRing(RichRingR* ring);
127 static float getCollectedPhotoElectrons();
128 static Double_t GetDistanceToTileBorder(Double_t x, Double_t y);
129 static int GetTileNumber(Double_t x,Double_t y);
130 
131 
132 myRichRingPlusFiller::myRichRingPlusFiller(){
133 #ifdef PDEBUG
134  printf("In myRichRingPlusFiller::myRichRingPlusFiller\n");
135 #endif
136  PRINTDEBUG;
137  init();
138  PRINTDEBUG;
139 }
140 
141 myRichRingPlusFiller::~myRichRingPlusFiller(){
142 #ifdef PDEBUG
143  printf("In myRichRingPlusFiller::~myRichRingPlusFiller\n");
144 #endif
145  PRINTDEBUG;
146  Clear();
147  PRINTDEBUG;
148  return;
149 }
150 
151 void myRichRingPlusFiller::Clear(Option_t* option){
152 #ifdef PDEBUG
153  printf("In myRichRingPlusFiller::Clear\n");
154 #endif
155  PRINTDEBUG;
156  ring=0;
157  Status = 999999;
158  PRINTDEBUG;
159  return;
160 }
161 
162 void myRichRingPlusFiller::init(){
163 #ifdef PDEBUG
164  printf("In myRichRingPlusFiller::init\n");
165 #endif
166  PRINTDEBUG;
167  Clear();
168  PRINTDEBUG;
169  return;
170 }
171 
172 void myRichRingPlusFiller::Fill(RichRingR* _ring, ParticleR* pp){
173 #ifdef PDEBUG
174  printf("In myRichRingPlusFiller::Fill\n");
175 #endif
176 
177  ring = _ring;
178 
179  if (ring && pp) {
180 
181  Status=ring->Status;
182 
183  if (pp->RichParticles<=1) JustOneRichParticle=true;//No more than 1 rich crossing particle
184 
185  if (ring->pTrTrack() && ring->pTrTrack()==pp->pTrTrack()) TheCorrectTrack=true;//è più della richiesta che fà carlitos nell'espempio in gbatch..
186 
187  if (ring->IsClean()) _IsClean=true;//Ring not contaminated by crossing particles
188  // if ((ringstatus&1)==0) {//Ring not contaminated by crossing particles
189 
190  Double_t r2=ring->TrRadPos[0]*ring->TrRadPos[0]+ring->TrRadPos[1]*ring->TrRadPos[1];
191 
192  if (r2<=60*60) FiducialVolumeLoose=true;//Crossing point within fiducial volume
193 
194  Double_t dist=GetDistanceToTileBorder(ring->TrRadPos[0], ring->TrRadPos[1]);
195 
196  if (dist>0.5) FiducialVolumeTight=true;//provare anche dist>1.0, Crossing point not close to tile border
197 
198  if (ring->Used>3) EnoughUsedHitsInRing=true;//Enough hits in ring
199 
200  if (ring->UsedM==0) NoMirroredPhotons=true;//nessun fotone è stato riflesso (provare ma solo se proprio necessario)
201 
202  if (ring->NpExp>2) EnoughExpetedHitsInRing=true;//Enough expected hits in ring
203 
204  if (ring->Prob>0.1) GoodRingProbability=true;//Ring hits distribution close to expected one
205 
206  if (ring->RingWidth(false)<3) RingWidthSmallEnough=true;//Ring Width small enough
207  // if (betaWidth(ring)<2.5) {//Ring Width small enough
208 
209  if (getPMTs()>4) EnoghPMTsInEvent=true;
210 
211  if(getPMTsForRing(ring)>2) EnoughPMTsInEventForRing=true;
212 
213  if(ring->getPhotoElectrons()>0.5*getCollectedPhotoElectrons()) EnoughPhotoElectronsWRTCollected=true;
214 
215  if (pp->pRichRingB()) {
216  if(fabs(pp->RichBetasDiscrepancy())<1e-3) BetaDiscrepancy=true;
217  // if (((ring->Beta-ringb->Beta)/(ring->Beta+ringb->Beta))>0.001) discr=true;
218  }
219 
220  PRINTDEBUG;
221 
222  Beta = ring->Beta;
223  Theta = ring->Theta;
224  Used = ring->Used;
225  UsedM = ring->UsedM;
226  Prob = ring->getProb();
227  copy_n( ring->AMSTrPars, 5,AMSTrPars);
228  copy_n( ring->TrRadPos , 3,TrRadPos);
229  copy_n( ring->TrPMTPos , 3,TrPMTPos);
230  BetaConsistency = ring->getBetaConsistency();
231  _betaCorrection = ring->betaCorrection();
232  PhotoElectrons = ring->getPhotoElectrons();
233  ExpectedPhotoElectrons = ring->getExpectedPhotoElectrons();
234  PMTs = ring->getPMTs();
235  PMTChargeConsistency = ring->getPMTChargeConsistency();
236  _PmtCorrectionsFailed = ring->PmtCorrectionsFailed();
237  TileIndex = ring->getTileIndex();
238  CollectedPhotoElectrons = RichHitR::getCollectedPhotoElectrons();
239  UsedHits = ring->getUsedHits(true);
240 
241  ChargeCorrections = ring->buildChargeCorrections();
242 
243  _RingWidth = ring->RingWidth();
244  BetaExpectedResolution = ring->getBetaExpectedResolution();
245  BetaExpectedRms = ring->getBetaExpectedRms();
246  ChargeExpectedResolution= ring->getChargeExpectedResolution();
247  ChargeExpectedRms = ring->getChargeExpectedRms();
248  NpExp = ring->NpExp;
249  NpCol = ring->NpCol;
250 
251  PRINTDEBUG;
252 
253  if(pp->pTrTrack()){
254  Double_t mindist[2], empoint[3], emdir[2];
255  AMSPoint emp;
256  AMSDir emd;
257  AMSPoint dist_p;
258  PRINTDEBUG;
259  kind = richrec( pp->pTrTrack(), emp, emd, dist_p, nphoel, npmthits, npmts, betam);
260  PRINTDEBUG;
261  empoint[0]=emp.x();
262  empoint[1]=emp.y();
263  empoint[2]=emp.z();
264  emdir[0]=emd.gettheta();
265  emdir[1]=emd.getphi();
266  mindist[0]=dist_p.x();
267  mindist[1]=dist_p.y();
268  }
269 
270  PRINTDEBUG;
271 
272  }
273 
274  PRINTDEBUG;
275 
276  return;
277 }
278 
279 Double_t betaWidth(RichRingR* ring){
280  Double_t weight=0;
281  Double_t sum=0;
282  for(int i=0;i<10;i++){
283  // cout<<" "<<i<<"=="<<ring->UsedWindow[i];
284  Double_t w=ring->UsedWindow[i]-(i>0?ring->UsedWindow[i-1]:0);
285  weight+=w;
286  sum+=(0.5+i)*(0.5+i)*w;
287  }
288  // cout<<" gives "<<sqrt(2*sum/weight)<<endl;
289  return sqrt(2*sum/weight);
290 }
291 
292 int GetTileNumber(Double_t x, Double_t y){
293  int nx=int(x/tile_width+5.5);
294  int ny=int(y/tile_width+5.5);
295  int tile=ny*grid_side_length+nx;
296  if(tile<0 || tile>n_tiles){
297  cout<<"GetTileNumber returns a nonsensical value: "<<endl
298  <<" x="<<x<<" y="<<y<<" nx= "<<nx<<" ny="<<ny<<" tile="<<tile<<endl;
299  return -1;
300  }
301  return tile;
302 }
303 
304 Double_t GetDistanceToTileBorder(Double_t x, Double_t y){
305  Double_t dx=fabs(fmod(x+tile_width/2,tile_width))-tile_width/2;
306  Double_t dy=fabs(fmod(y+tile_width/2,tile_width))-tile_width/2;
307  return tile_width/2-(fabs(dx+dy)/2+fabs(dx-dy)/2);
308 }
309 
310 int getPMTs(bool countCrossed){
311  AMSEventR* event = AMSEventR::Head();
312  if(!event) return 0;
313  bool counted[680];
314  for(int i=0;i<680;counted[i++]=false);
315  int counter=0;
316  for(int i=0;i<event->nRichHit();i++){
317  RichHitR *hit=event->pRichHit(i);
318  if(!hit) continue;
319  if(!countCrossed && hit->IsCrossed()) continue;
320  int pmt=hit->Channel/16;
321  if(counted[pmt]) continue;
322  counted[pmt]=true;
323  counter++;
324  }
325  return counter;
326 }
327 
328 int getPMTsForRing(RichRingR* ring){
329  AMSEventR* event = AMSEventR::Head();
330  if(!event) return 0;
331  bool counted[680];
332  for(int i=0;i<680;counted[i++]=false);
333  int counter=0;
334  for(int i=0;i<ring->Used;i++){
335  RichHitR *hit=event->pRichHit(ring->iRichHit(i));
336  if(!hit) continue;
337  if(hit->IsCrossed()) continue;
338  int pmt=hit->Channel/16;
339  if(counted[pmt]) continue;
340  counted[pmt]=true;
341  counter++;
342  }
343  return counter;
344 }
345 
346 float getCollectedPhotoElectrons(){
347  AMSEventR* event = AMSEventR::Head();
348  if(!event) return 0;
349  float counter=0;
350  for(int i=0;i<event->nRichHit();i++){
351  RichHitR *hit=event->pRichHit(i);
352  if(!hit) continue;
353  if(hit->IsCrossed()) continue;
354  counter+=hit->Npe;
355  }
356  return counter;
357 }
358 
359 #endif //#ifdef _WITHGBATCH_
Generated on Tue Dec 29 2015 23:02:50 for AMSDST by   doxygen 1.8.3