BOSS 6.6.4.p01
BESIII Offline Software System
Loading...
Searching...
No Matches
TofShower.cxx
Go to the documentation of this file.
1#include "GaudiKernel/MsgStream.h"
2#include "GaudiKernel/Bootstrap.h"
3#include "GaudiKernel/ISvcLocator.h"
4#include "GaudiKernel/PropertyMgr.h"
5#include "GaudiKernel/IJobOptionsSvc.h"
6#include "GaudiKernel/INTupleSvc.h"
7#include "TofEnergyRec/TofShower.h"
8#include "Identifier/Identifier.h"
9#include "Identifier/TofID.h"
10#include "RawDataProviderSvc/TofData.h"
11#include "TofCaliSvc/ITofCaliSvc.h"
12#include "TofQCorrSvc/ITofQCorrSvc.h"
13#include "TofQElecSvc/ITofQElecSvc.h"
14#include "DstEvent/TofHitStatus.h"
15#include "globals.hh"
16#include <G4String.hh>
17#include <iostream>
18#include <fstream>
19#include <math.h>
20#include <string>
21#include "TofSim/BesTofDigitizer.hh"
22#include "EvTimeEvent/RecEsTime.h"
23#include "RawEvent/RawDataUtil.h"
24
25TofShower::TofShower():m_output(false),m_isData(true) // m_output wird auf false und m_isData auf true gesetzt
26{
27 IJobOptionsSvc* jobSvc;
28 Gaudi::svcLocator()->service("JobOptionsSvc", jobSvc);
29 jobSvc->setMyProperties("TofShower", &m_propMgr);
30
31}
32
33void TofShower::BookNtuple(NTuple::Tuple*& tuple, NTuple::Tuple*& tuple1, NTuple::Tuple*& tuple2)
34{
35 m_output = true;
36 m_tuple = tuple;
37 if(!m_tuple) {
38 std::cerr << "Invalid ntuple in TofEnergyRec!" << std::endl;
39 } else {
40 m_tuple->addItem ("part", m_part);
41 m_tuple->addItem ("layer", m_layer);
42 m_tuple->addItem ("im", m_im);
43 m_tuple->addItem ("end", m_end);
44 m_tuple->addItem ("zpos", m_zpos);
45 m_tuple->addItem ("adc1", m_adc1);
46 m_tuple->addItem ("adc2", m_adc2);
47 m_tuple->addItem ("tdc1", m_tdc1);
48 m_tuple->addItem ("tdc2", m_tdc2);
49 m_tuple->addItem ("energy", m_energy);
50 }
51
52 m_tuple1 = tuple1;
53 if(!m_tuple1) {
54 std::cerr << "Invalid ntuple1 in TofEnergyRec!" << std::endl;
55 } else {
56 m_tuple1->addItem ("part", m_shower_part);
57 m_tuple1->addItem ("layer", m_shower_layer);
58 m_tuple1->addItem ("im", m_shower_im);
59 m_tuple1->addItem ("zpos", m_shower_zpos);
60 m_tuple1->addItem ("energy", m_shower_energy);
61 }
62
63 m_tuple2 = tuple2;
64 if(!m_tuple2) {
65 std::cerr << "Invalid ntuple2 in TofEnergyRec!" << std::endl;
66 } else {
67 m_tuple2->addItem ("dist", m_seed_dist);
68 }
69}
70
71void TofShower::energyCalib(vector<TofData*>& tofDataVec, RecTofTrackCol* recTofTrackCol)
72{
73 //Get TOF Calibtration Service
74 ISvcLocator* svcLocator = Gaudi::svcLocator();
75 ITofCaliSvc* tofCaliSvc;
76 StatusCode sc = svcLocator->service("TofCaliSvc", tofCaliSvc);
77 if (sc != StatusCode::SUCCESS) {
78 cout << "TofEnergyRec Get Calibration Service Failed !! " << endl;
79 }
80
81 ITofQCorrSvc* tofQCorrSvc;
82 sc = svcLocator->service("TofQCorrSvc", tofQCorrSvc);
83 if (sc != StatusCode::SUCCESS) {
84 cout << "TofEnergyRec Get QCorr Service Failed !! " << endl;
85 }
86
87 ITofQElecSvc* tofQElecSvc;
88 sc = svcLocator->service("TofQElecSvc", tofQElecSvc);
89 if (sc != StatusCode::SUCCESS) {
90 cout << "TofEnergyRec Get QElec Service Failed !! " << endl;
91 }
92
93 vector<TofData*>::iterator it;
94 for(it=tofDataVec.begin();
95 it!=tofDataVec.end();
96 it++) {
97
98 Identifier id((*it)->identify());
99 int barrel_ec = TofID::barrel_ec(id);
100 int layer = TofID::layer(id);
101 int im = TofID::phi_module(id);
102 int end = TofID::end(id);
103 bool is_mrpc = TofID::is_mymrpc(id);
104
105
106 if(m_output) {
107 m_part = barrel_ec;
108 m_layer = layer;
109 m_im = im;
110 m_end = end;
111 }
112
113 if((*it)->barrel()) {
114 TofData* bTofData = (*it);
115 bTofData->setZpos(99.);
116 bTofData->setEnergy(0.);
117 if(bTofData->tdc1()<=0||bTofData->tdc1()>8000||bTofData->tdc2()<=0||bTofData->tdc2()>8000) continue;
118
119 double adc1,adc2,tdc1,tdc2;
120 tdc1 = bTofData->tdc1();
121 tdc2 = bTofData->tdc2();
122 adc1 = bTofData->adc1();
123 adc2 = bTofData->adc2();
124
125 //from data CalibSvc
126 double zpos = tofCaliSvc->ZTDC( tdc1, tdc2, bTofData->tofId() );
127 if(fabs(zpos)>115) continue;
128 double tofq = tofCaliSvc->BPh( adc1, adc2, zpos, bTofData->tofId());
129 if(tofq<100||tofq>10000) continue;
130 //double energy = q*0.0036;
131 double energy = tofq*m_calibConst; //new calibration result in 2009.9.27
132 //cout<< "TofShower Barrel energy " << energy << endl;
133 zpos /= 100.; //cm->m
134
135 bTofData->setZpos(zpos);
136 bTofData->setEnergy(energy);
137
138 if(m_output) {
139 m_part = barrel_ec;
140 m_layer = layer;
141 m_im = im;
142 m_end = end;
143 m_adc1 = bTofData->adc1();
144 m_adc2 = bTofData->adc2();
145 m_tdc1 = bTofData->tdc1();
146 m_tdc2 = bTofData->tdc2();
147 m_zpos = zpos;
148 m_energy = energy;
149 m_tuple->write();
150 }
151
152 } else {
153
154 if(!is_mrpc)
155 {
156 //cout<<"endcap"<<endl;
157 //ETofData* eTofData = dynamic_cast<ETofData*>(*it);
158 //TofData* bTofData = (*it);
159 //double energy = 2*eTofData->adcChannel()/140;
160 //eTofData->setEnergy(energy);
161 }
162 else if(is_mrpc)
163 {
164 //cout<< "TofShower : Energy reconstruction for the MRPC detector "<<endl;
165 TofData* eTofData = (*it);
166 double charge_ns = eTofData->adc(); // time over threshold in ns
167 //cout<< " partID : "<< barrel_ec << endl;
168 //cout<< " module : "<<im/25 << endl;
169 //cout<< " deposit charge adccha : " << eTofData->adcChannel() <<endl;
170 //cout<< " deposit charge adc : " << eTofData->adc() <<endl;
171 //cout<< "_____________________________________" << endl;
172
173 //The next step is to determine the energy deposit in the module. We have three different cases: Neutral tracks, Charged tracks with neighbour, charged track without neighbour
174 bool neutral =true;
175 bool neighborhood = false; //Are more than two cells fired?
176
177 int partId_all = barrel_ec;
178 int module_all = im/25;
179 int strip_all = im%25;
180
181 int partId_charged =-9999;
182 int module_charged =-9999;
183 int strip_charged =-9999;
184
185
186
187
188 /*
189 cout<< "TofShower : All modules "<<endl;
190 cout << "partid_all " << partId_all<<endl;
191 cout << "module_all " << module_all<<endl;
192 cout << "strip_all " << strip_all << endl;
193 cout<< "_____________________________________" << endl;
194 */
195
196
197 RecTofTrackCol::iterator rec_charged;
198
199
200 for( rec_charged = recTofTrackCol->begin(); rec_charged != recTofTrackCol->end();rec_charged++)
201 {
202
203
204
205
206
207 TofHitStatus *status = new TofHitStatus;
208 status->setStatus((*rec_charged)->status());
209 if((status->is_barrel())) continue; //Consider only endcap tracks
210 if((*rec_charged)->tofID()<=125) continue;
211
212 if(status->is_east())
213 {
214 if((*rec_charged)->tofID()>=600) {partId_charged=3; module_charged = ((*rec_charged)->tofID()-574)/25; strip_charged =((*rec_charged)->tofID()-574)%25; }
215 else {partId_charged=4; module_charged = ((*rec_charged)->tofID()-100)/25; strip_charged =((*rec_charged)->tofID()-100)%25; }
216 }
217 else
218 {
219 if((*rec_charged)->tofID()>=600) {partId_charged=6; module_charged = ((*rec_charged)->tofID()-574)/25; strip_charged =((*rec_charged)->tofID()-574)%25; }
220 else {partId_charged=5;module_charged = ((*rec_charged)->tofID()-100)/25; strip_charged =((*rec_charged)->tofID()-100)%25;}
221 }
222
223 delete status;
224
225 /*
226 cout<< " TofShower : Charged modules "<<endl;
227 cout << " partid_ch " << partId_charged<<endl;
228 cout << " module_ch " << module_charged<<endl;
229 cout << " strip_ch " << strip_charged << endl;
230 */
231
232 //Determine wether a track is neutral or not:
233 if( (partId_all==3 || partId_all==4) && (partId_charged==3 || partId_charged==4 ) )
234 {
235 if( partId_all==partId_charged && module_charged == module_all && abs(strip_charged-strip_all)<3) neutral =false;
236 if( partId_all==partId_charged && module_charged == module_all)
237 {
238 if(strip_all%2==0 && strip_charged%2!=0 && strip_charged<strip_all && abs(strip_charged-strip_all)<4) neutral =false;
239 if(strip_all%2!=0 && strip_charged%2==0 && strip_charged>strip_all && abs(strip_charged-strip_all)<4) neutral =false;
240 }
241
242
243 if( partId_all!=partId_charged && abs(module_charged-module_all)<2)
244 {
245 if(strip_all%2==0 && strip_charged%2!=0 && abs(strip_charged-strip_all)<3 ) neutral =false;
246 if(strip_all%2==0 && strip_charged%2!=0 && strip_charged<strip_all && abs(strip_charged-strip_all)<4) neutral =false;
247
248 if(strip_all%2!=0 && strip_charged%2==0 && abs(strip_charged-strip_all)<3 ) neutral =false;
249 if(strip_all%2!=0 && strip_charged%2==0 && strip_charged>strip_all && abs(strip_charged-strip_all)<4) neutral =false;
250 }
251
252 }
253 else if ( (partId_all==5 || partId_all==6) && (partId_charged==5 || partId_charged==6 ) )
254 {
255 if( partId_all==partId_charged && module_charged ==module_all && abs(strip_charged-strip_all)<3) neutral =false;
256
257 if( partId_all==partId_charged && module_charged == module_all)
258 {
259 if(strip_all%2==0 && strip_charged%2!=0 && strip_charged<strip_all && abs(strip_charged-strip_all)<4) neutral =false;
260 if(strip_all%2!=0 && strip_charged%2==0 && strip_charged>strip_all && abs(strip_charged-strip_all)<4) neutral =false;
261 }
262
263
264 if( partId_all!=partId_charged && abs(module_charged-module_all)<2)
265 {
266 if(strip_all%2==0 && strip_charged%2!=0 && abs(strip_charged-strip_all)<3 ) neutral =false;
267 if(strip_all%2==0 && strip_charged%2!=0 && strip_charged<strip_all && abs(strip_charged-strip_all)<4) neutral =false;
268
269 if(strip_all%2!=0 && strip_charged%2==0 && abs(strip_charged-strip_all)<3 ) neutral =false;
270 if(strip_all%2!=0 && strip_charged%2==0 && strip_charged>strip_all && abs(strip_charged-strip_all)<4) neutral =false;
271 }
272 }
273
274 /*
275 cout<< "TofShower : Energy reconstruction for the MRPC detector "<<endl;
276 cout << "->tofID() " <<(*rec_charged)->tofID()<<endl;
277 cout << "->status() " <<(*rec_charged)->status()<<endl;
278 cout << "partid_charged " << partId_charged<<endl;
279 cout << "module_charged " << module_charged<<endl;
280 cout << "strip_charged " << strip_charged << endl;
281 cout << "neutral ? " << neutral << endl;
282 cout<< "_____________________________________" << endl;
283 */
284 }//Close for rec charegd Track
285
286 //Yet we do know wether the track is neutral or not
287
288 //Next step: Determine wether the track has neighborhood (Have two strips in a neighborhood fired?)
289 vector<TofData*>::iterator it_temp;
290 for(it_temp=tofDataVec.begin();it_temp!=tofDataVec.end();it_temp++)
291 {
292
293 Identifier id_temp((*it_temp)->identify());
294 int partId_temp = TofID::barrel_ec(id_temp);
295 int im_temp = TofID::phi_module(id_temp);
296 int module_temp = im_temp/25;
297 int strip_temp = im_temp%25;
298 bool is_mrpc_temp = TofID::is_mymrpc(id_temp);
299
300 if( partId_temp==1 ) continue;
301 if(!is_mrpc_temp) continue;
302
303 if( (partId_all==3 || partId_all==4) && (partId_temp==3 || partId_temp==4 ) )
304 {
305 if( partId_all==partId_temp && module_temp == module_all && abs(strip_temp-strip_all)!=0 && abs(strip_temp-strip_all)<3) neighborhood =true;
306 if( partId_all==partId_temp && module_temp == module_all)
307 {
308 if(strip_all%2==0 && strip_temp%2!=0 && strip_temp<strip_all && abs(strip_temp-strip_all)<4 && abs(strip_temp-strip_all)!=0) neighborhood =true;
309 if(strip_all%2!=0 && strip_temp%2==0 && strip_temp>strip_all && abs(strip_temp-strip_all)<4 && abs(strip_temp-strip_all)!=0) neighborhood =true;
310 }
311
312 if( partId_all!=partId_temp && abs(module_temp-module_all)<2)
313 {
314 if(strip_all%2==0 && strip_temp%2!=0 && abs(strip_temp-strip_all)<3 ) neighborhood =true;
315 if(strip_all%2==0 && strip_temp%2!=0 && strip_temp<strip_all && abs(strip_temp-strip_all)<4) neighborhood =true;
316
317 if(strip_all%2!=0 && strip_temp%2==0 && abs(strip_temp-strip_all)<3 ) neighborhood =true;
318 if(strip_all%2!=0 && strip_temp%2==0 && strip_temp>strip_all && abs(strip_temp-strip_all)<4) neighborhood =true;
319 }
320
321 }//close if partid 3 & 4
322 else if( (partId_all==5 || partId_all==6) && (partId_temp==5 || partId_temp==6 ) )
323 {
324 if( partId_all==partId_temp && module_temp == module_all && abs(strip_temp-strip_all)!=0 && abs(strip_temp-strip_all)<3) neighborhood =true;
325 if( partId_all==partId_temp && module_temp == module_all)
326 {
327 if(strip_all%2==0 && strip_temp%2!=0 && strip_temp<strip_all && abs(strip_temp-strip_all)<4 && abs(strip_temp-strip_all)!=0) neighborhood =true;
328 if(strip_all%2!=0 && strip_temp%2==0 && strip_temp>strip_all && abs(strip_temp-strip_all)<4 && abs(strip_temp-strip_all)!=0) neighborhood =true;
329 }
330
331 if( partId_all!=partId_temp && abs(module_temp-module_all)<2)
332 {
333 if(strip_all%2==0 && strip_temp%2!=0 && abs(strip_temp-strip_all)<3 ) neighborhood =true;
334 if(strip_all%2==0 && strip_temp%2!=0 && strip_temp<strip_all && abs(strip_temp-strip_all)<4) neighborhood =true;
335
336 if(strip_all%2!=0 && strip_temp%2==0 && abs(strip_temp-strip_all)<3 ) neighborhood =true;
337 if(strip_all%2!=0 && strip_temp%2==0 && strip_temp>strip_all && abs(strip_temp-strip_all)<4) neighborhood =true;
338 }
339
340 }//close if partid 5 & 6
341
342 }
343
344
345 double energy =0; // in MeV
346
347 if(neutral==true)
348 {
349 // The values are from photons simulated with geant4!
350
351 if(charge_ns<15.0716) energy = 326.35-10.6284*15.0716-1.4169*15.0716*15.0716-0.0514292*15.0716*15.0716*15.0716+0.00645421*15.0716*15.0716*15.0716*15.0716;//in MeV: Below this value we assume a constant energyloss!
352 else if(charge_ns>15.9131) energy = -7.91205+0.814588*charge_ns; //Over this value we expect the energyloss to be the same as for charged particles!
353 else energy = 326.35-10.6284*charge_ns-1.4169*charge_ns*charge_ns-0.0514292*charge_ns*charge_ns*charge_ns+0.00645421*charge_ns*charge_ns*charge_ns*charge_ns;// in MeV
354 }
355 else //charged tracks, (combined electrons/positrons in geant4)
356 {
357 if(neighborhood==false) energy = -7.91205+0.814588*charge_ns;// in MeV
358 else
359 {
360 if(15.16>= charge_ns) energy= -72550.7+14230.9*15.16-930.499*15.16*15.16+20.2818*15.16*15.16*15.16;//Minimal loss from fit!
361 if(15.16< charge_ns &&charge_ns<15.7633) energy= -72550.7+14230.9*charge_ns-930.499*charge_ns*charge_ns+20.2818*charge_ns*charge_ns*charge_ns;
362 else energy = -7.91205+0.814588*charge_ns; //Else we will use the not neighborhood distribution as our fit will begin to fall
363 }
364 }//Close else neutral
365
366
367
368 eTofData->setEnergy(energy);//set the energy loss in MeV
369
370 /*
371 cout << "-----------------" << endl;
372 cout<< "TofShower : Set the energy for the individual TofTracks" << endl;
373 cout << "partid_all " << partId_all<<endl;
374 cout << "module_all " << module_all<<endl;
375 cout << "strip_all " << strip_all << endl;
376 cout << "Pulselenght " << charge_ns <<endl;
377 cout << "neutral ? " << neutral << endl;
378 cout << "neighborhood ? " << neighborhood<< endl;
379 cout << "Energy [MeV] " << energy<< endl;
380 cout<< "_____________________________________" << endl;
381 */
382
383
384
385 }//close if(is_mrpc)
386
387 }//close else (endcap)
388 }
389}
390
391vector<Identifier> TofShower::getNeighbors(const Identifier& id)
392{
393 vector<int> NeighborVec;
394 vector<Identifier> NeighborIdVec;
395 NeighborVec.clear();
396 NeighborIdVec.clear();
397
398 int barrel_ec = TofID::barrel_ec(id);
399 int layer = TofID::layer(id);
400 int im = TofID::phi_module(id);
401 int end = TofID::end(id);
402
403 if(barrel_ec==1) { //barrel
404 int num = im+layer*88;
405 if(num<88) { //layer1
406 if(num==0) {
407 NeighborVec.push_back(1);
408 NeighborVec.push_back(87);
409 NeighborVec.push_back(88);
410 NeighborVec.push_back(89);
411 } else if(num==87) {
412 NeighborVec.push_back(0);
413 NeighborVec.push_back(86);
414 NeighborVec.push_back(88);
415 NeighborVec.push_back(175);
416 } else {
417 NeighborVec.push_back(num+1);
418 NeighborVec.push_back(num-1);
419 NeighborVec.push_back(num+88);
420 NeighborVec.push_back(num+88+1);
421 }
422 } else {
423 if(num==88) {
424 NeighborVec.push_back(89);
425 NeighborVec.push_back(175);
426 NeighborVec.push_back(0);
427 NeighborVec.push_back(87);
428 } else if(num==175) {
429 NeighborVec.push_back(88);
430 NeighborVec.push_back(174);
431 NeighborVec.push_back(86);
432 NeighborVec.push_back(87);
433 } else {
434 NeighborVec.push_back(num+1);
435 NeighborVec.push_back(num-1);
436 NeighborVec.push_back(num-88);
437 NeighborVec.push_back(num-88-1);
438 }
439 }
440
441 int size=NeighborVec.size();
442 for(int i=0;i<size;i++) {
443 layer = NeighborVec[i]/88;
444 im = NeighborVec[i]%88;
445 NeighborIdVec.push_back(TofID::cell_id(barrel_ec,layer,im,end));
446 }
447 }
448
449 else { //endcap
450 if(im==0) {
451 NeighborVec.push_back(1);
452 NeighborVec.push_back(47);
453 } else if(im==47) {
454 NeighborVec.push_back(0);
455 NeighborVec.push_back(46);
456 } else {
457 NeighborVec.push_back(im-1);
458 NeighborVec.push_back(im+1);
459 }
460
461 int size=NeighborVec.size();
462 for(int i=0;i<size;i++) {
463 im = NeighborVec[i];
464 NeighborIdVec.push_back(TofID::cell_id(barrel_ec,layer,im,end));
465 }
466 }
467
468 return NeighborIdVec;
469}
470
471vector<Identifier> TofShower::getNextNeighbors(const Identifier& id)
472{
473 vector<Identifier> NeighborVec,tmpNeighborVec,tmpNextNeighborVec;
474 vector<Identifier>::iterator ci_NV,ci_tmpNV,ci_tmpNNV;
475 NeighborVec=getNeighbors(id);
476 tmpNeighborVec=getNeighbors(id);
477 bool flag=false; //whether NeighborVec already includes this crystal
478 bool flagNeighbor=false; //whether this crystal belongs to NeighborVec
479
480 //------------------------------------------------------------------
481 for(ci_tmpNV=tmpNeighborVec.begin();
482 ci_tmpNV!=tmpNeighborVec.end();
483 ci_tmpNV++){
484 tmpNextNeighborVec=getNeighbors(*ci_tmpNV);
485 //================================================================
486 for(ci_tmpNNV=tmpNextNeighborVec.begin();
487 ci_tmpNNV!=tmpNextNeighborVec.end();
488 ci_tmpNNV++){
489
490 for(ci_NV=NeighborVec.begin();
491 ci_NV!=NeighborVec.end();
492 ci_NV++){
493 if(*ci_tmpNNV==*ci_NV){ //this crystal is already included
494 flag=true;
495 break;
496 }
497 }
498
499 if(!flag){ //find a new crystal
500 //for(ci_tmpNV1=tmpNeighborVec.begin();
501 // ci_tmpNV1!=tmpNeighborVec.end();
502 // ci_tmpNV1++){
503 // if(*ci_tmpNNV==*ci_tmpNV1){ //this crystal belongs to NeighborVec
504 // flagNeighbor=true;
505 // break;
506 // }
507 //}
508
509 if(*ci_tmpNNV==id)
510 flagNeighbor=true;
511
512 if(!flagNeighbor)
513 NeighborVec.push_back(*ci_tmpNNV);
514 else
515 flagNeighbor=false;
516 }
517 else
518 flag=false;
519 }
520 //================================================================
521 }
522 //------------------------------------------------------------------
523
524 return NeighborVec;
525}
526
527
528vector<Identifier> TofShower::getMRPC_neighbours(const Identifier& id)
529{
530 //This function returns the following neigbours (x=base)
531 // - - -
532 // - - -
533 // - x -
534 // - - -
535 // - - -
536
537
538 vector<Identifier> NeighborVec;
539
540 int partid_it = TofID::barrel_ec(id);
541 int tofid_it = TofID::phi_module(id);
542
543
544 //int module= BesTofDigitizerEcV4::Get_module_mrpc_from_unique_identifier(partid_it);
545
546 int tofid_l = get_mrpc_stripid_neighbour("l",tofid_it,partid_it);
547 int partid_l = get_mrpc_partid_neighbour("l",tofid_it,partid_it);
548 int tofid_r = get_mrpc_stripid_neighbour("r",tofid_it,partid_it);
549 int partid_r = get_mrpc_partid_neighbour("r",tofid_it,partid_it);
550 int tofid_u = get_mrpc_stripid_neighbour("u",tofid_it,partid_it);
551 int partid_u =get_mrpc_partid_neighbour("u",tofid_it,partid_it);
552 int tofid_b = get_mrpc_stripid_neighbour("b",tofid_it,partid_it);
553 int partid_b = get_mrpc_partid_neighbour("b",tofid_it,partid_it);
554
555 int tofid_lu = get_mrpc_stripid_neighbour("u",tofid_l,partid_l);
556 int partid_lu = get_mrpc_partid_neighbour("u",tofid_l,partid_l);
557 int tofid_lb = get_mrpc_stripid_neighbour("b",tofid_l,partid_l);
558 int partid_lb = get_mrpc_partid_neighbour("b",tofid_l,partid_l);
559 int tofid_ru = get_mrpc_stripid_neighbour("u",tofid_r,partid_r);
560 int partid_ru = get_mrpc_partid_neighbour("u",tofid_r,partid_r);
561 int tofid_rb = get_mrpc_stripid_neighbour("b",tofid_r,partid_r);
562 int partid_rb = get_mrpc_partid_neighbour("b",tofid_r,partid_r);
563
564 int tofid_luu = get_mrpc_stripid_neighbour("u",tofid_lu,partid_lu);
565 int partid_luu =get_mrpc_partid_neighbour("u",tofid_lu,partid_lu);
566 int tofid_lbb = get_mrpc_stripid_neighbour("b",tofid_lb,partid_lb);
567 int partid_lbb = get_mrpc_partid_neighbour("b",tofid_lb,partid_lb);
568 int tofid_ruu = get_mrpc_stripid_neighbour("u",tofid_ru,partid_ru);
569 int partid_ruu = get_mrpc_partid_neighbour("u",tofid_ru,partid_ru);
570 int tofid_rbb = get_mrpc_stripid_neighbour("b",tofid_rb,partid_rb);
571 int partid_rbb = get_mrpc_partid_neighbour("b",tofid_rb,partid_rb);
572 int tofid_uu = get_mrpc_stripid_neighbour("u",tofid_u,partid_u);
573 int partid_uu = get_mrpc_partid_neighbour("u",tofid_u,partid_u);
574 int tofid_bb = get_mrpc_stripid_neighbour("b",tofid_b,partid_b);
575 int partid_bb = get_mrpc_partid_neighbour("b",tofid_b,partid_b);
576
577
578 /*
579 cout << "TofShower::getMRPC_neighbours: " << endl;
580 cout << " IT " << partid_it << " | " << tofid_it <<endl;
581 cout << " l " << partid_l<< " | " << tofid_l <<endl;
582 cout << " r " << partid_r<< " | " << tofid_r <<endl;
583 cout << " u " << partid_u<< " | " << tofid_u <<endl;
584 cout << " b " << partid_b<< " | " << tofid_b <<endl;
585
586 cout << " lu " << partid_lu<< " | " << tofid_lu <<endl;
587 cout << " lb " << partid_lb<< " | " << tofid_lb <<endl;
588 cout << " ru " << partid_ru<< " | " << tofid_ru <<endl;
589 cout << " rb " << partid_rb<< " | " << tofid_rb <<endl;
590
591 cout << " luu " << partid_luu<< " | " << tofid_luu <<endl;
592 cout << " lbb " << partid_lbb<< " | " << tofid_lbb <<endl;
593 cout << " ruu " << partid_ruu<< " | " << tofid_ruu <<endl;
594 cout << " rbb " << partid_rbb<< " | " << tofid_rbb <<endl;
595 cout << " uu " << partid_uu<< " | " << tofid_uu <<endl;
596 cout << " bb " << partid_bb<< " | " << tofid_bb <<endl;
597 */
598
599
600
601
602 //if(tofid_it%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_it,tofid_it));
603 if(tofid_l%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_l,tofid_l));
604 if(tofid_r%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_r,tofid_r));
605 if(tofid_u%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_u,tofid_u));
606 if(tofid_b%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_b,tofid_b));
607
608 if(tofid_lu%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_lu,tofid_lu));
609 if(tofid_lb%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_lb,tofid_lb));
610 if(tofid_ru%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_ru,tofid_ru));
611 if(tofid_rb%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_rb,tofid_rb));
612
613 if(tofid_luu%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_luu,tofid_luu));
614 if(tofid_lbb%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_lbb,tofid_lbb));
615 if(tofid_ruu%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_ruu,tofid_ruu));
616 if(tofid_rbb%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_rbb,tofid_rbb));
617 if(tofid_uu%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_uu,tofid_uu));
618 if(tofid_bb%25!=0) NeighborVec.push_back(TofID::cell_id_mrpc(partid_bb,tofid_bb));
619
620
621 return NeighborVec;
622
623}//close getMRPC Neighbours
624
625
626void TofShower::findSeed(vector<TofData*>& tofDataVec)
627{
628 bool max=false;
629 m_seedVec.clear();
630
631 vector<TofData*>::iterator it;
632 for(it=tofDataVec.begin();
633 it!=tofDataVec.end();
634 it++) {
635 if((*it)->barrel()) { //barrel
636 TofData* bTofData = (*it);
637
638
639
640 //cout << "Identifier(bTofData->identify()) " << Identifier(bTofData->identify()) << endl;
641 //std::cout << "TofShower bTofData->energy() = " << bTofData->energy() << std::endl;
642 if(bTofData->energy()<5.) continue; //seed energy cut = 6MeV
643
644 max=true;
645 vector<Identifier> NeighborVec=getNextNeighbors(Identifier(bTofData->identify()));
646
647 //const Identifier & help = Identifier(bTofData->identify());
648 //cout << "TofShower::findSeed Barrel Track base "<<TofID::layer(help) << " " << TofID::phi_module(help) << endl;
649
650 vector<Identifier>::iterator iNeigh;
651 for(iNeigh=NeighborVec.begin();
652 iNeigh!=NeighborVec.end();
653 iNeigh++) {
654
655 //const Identifier & help2 = Identifier(*iNeigh);
656 //cout << "TofShower::findSeed Barrel Track neigh "<<TofID::layer(help2) << " " << TofID::phi_module(help2) << endl;
657
658 vector<TofData*>::iterator it2;
659 for(it2=tofDataVec.begin();
660 it2!=tofDataVec.end();
661 it2++) {
662 if((*it2)->identify()==*iNeigh) {
663 TofData* bTofData2 = (*it2);
664 if(bTofData2->energy()>bTofData->energy()) {
665 max=false;
666 }
667 break;
668 }
669 }
670
671 }
672 }
673 else { //endcap
674
675 //cout << "TofShower::findSeed Endcap Track" << endl;
676 TofData* eTofData = (*it);
677
678 Identifier id((*it)->identify());
679 bool is_mrpc = TofID::is_mymrpc(id);
680
681
682 if(!is_mrpc)
683 {
684 if(eTofData->energy()<5.) continue; //seed energy cut = 5MeV
685 max=true;
686 vector<Identifier> NeighborVec=getNextNeighbors(Identifier(eTofData->identify()));
687 vector<Identifier>::iterator iNeigh;
688 for(iNeigh=NeighborVec.begin(); iNeigh!=NeighborVec.end(); iNeigh++)
689 {
690
691 vector<TofData*>::iterator it2;
692 for(it2=tofDataVec.begin(); it2!=tofDataVec.end(); it2++)
693 {
694 if((*it2)->identify()==*iNeigh) {
695 TofData* eTofData2 = (*it2);
696 if(eTofData2->energy()>eTofData->energy()) {
697 max=false;
698 }
699 break;
700 }
701 }//close for
702 }//Close for
703 }//close if(!is_mrpc)
704 else if(is_mrpc)
705 {
706
707 //cout << "TofShower::findSeed MRPC Track "<< Identifier(eTofData->identify()) << endl;
708 if(eTofData->energy()<2.) continue; //Seed energy cut = 2 MeV
709 max=true;
710
711 vector<Identifier> NeighborVec=getMRPC_neighbours(Identifier(eTofData->identify()));
712
713 //cout << "TofShower::findSeed getMRPC_neighbours OK" << endl;
714 vector<Identifier>::iterator iNeigh;
715
716 for(iNeigh=NeighborVec.begin(); iNeigh!=NeighborVec.end(); iNeigh++)
717 {
718 //cout << "TofShower::findSeed Phimodule MRPC " << TofID::phi_module(*iNeigh) <<endl;
719 vector<TofData*>::iterator it2;
720 for(it2=tofDataVec.begin(); it2!=tofDataVec.end(); it2++)
721 {
722 if((*it2)->identify()==*iNeigh) {
723 TofData* eTofData2 = (*it2);
724 if(eTofData2->energy()>eTofData->energy()) {
725 max=false;
726 }
727 break;
728 }
729 }//close for
730 }//Close for
731 //cout << "TofShower::findSeed Both forloops done" << endl;
732
733 }//Close else if(is_mrpc)
734
735 } //close else endcap
736
737
738 if(max) {
739 m_seedVec.push_back(Identifier((*it)->identify()));
740 }
741
742 }//close loop over tofdata
743}//close find seed
744
745void TofShower::findShower(vector<TofData*>& tofDataVec, RecTofTrackCol* recTofTrackCol, double t0)
746{
747
748 energyCalib(tofDataVec, recTofTrackCol);
749 //cout << "TofShower::findShower energycalib OK" << endl;
750 findSeed(tofDataVec);
751 //cout << "TofShower::findShower findseed OK" << endl;
752
753 vector<Identifier>::iterator iSeed;
754
755 for(iSeed=m_seedVec.begin(); iSeed!=m_seedVec.end(); iSeed++)
756 {
757
758 bool is_mrpc = TofID::is_mymrpc(*iSeed);
759
760
761 if(!is_mrpc) //Old TofDetector + barrel
762 {
763 int barrel_ec = TofID::barrel_ec(*iSeed);
764 int layer = TofID::layer(*iSeed);
765 int im = TofID::phi_module(*iSeed);
766 im += layer * 88;
767
768 bool neutral=true;
769 //match with Tof charged track
770 int dphi=999;
771 RecTofTrackCol::iterator iTrack, iMatch;
772 for(iTrack=recTofTrackCol->begin(); iTrack!=recTofTrackCol->end(); iTrack++)
773 {
774 if(barrel_ec==1) {
775 dphi=abs(im-(*iTrack)->tofID());
776 dphi = dphi>=44 ? 88-dphi : dphi;
777 } else if(barrel_ec==2) {
778 dphi=abs(im-(*iTrack)->tofID()+48);
779 dphi = dphi>=24 ? 48-dphi : dphi;
780 } else {
781 dphi=abs(im-(*iTrack)->tofID());
782 dphi = dphi>=24 ? 48-dphi : dphi;
783 }
784 if(abs(dphi)<=2) {
785 iMatch = iTrack;
786 neutral=false;
787 break;
788 }
789 }
790
791 //energy sum of seed and its neighbors
792 //use avarage mean to calculation position
793 double zpos=0;
794 double energy=0;
795 double seedPos=0;
796
797
798 //cout << "Energhy =0 " << endl;
799 vector<TofData*>::iterator it;
800 for(it=tofDataVec.begin(); it!=tofDataVec.end(); it++)
801 {
802 if((*it)->identify()==*iSeed) {
803 //cout<<"iSeed="<<*iSeed<<endl;
804 if((*it)->barrel()) {
805 TofData* bTofData = (*it);
806 zpos+=bTofData->zpos()*bTofData->energy();
807 energy+=bTofData->energy();
808 seedPos=bTofData->zpos();
809
810 //cout << "Add energy barrel seed "<< TofID::layer(*iSeed) <<" " << TofID::phi_module(*iSeed) << " " << bTofData->energy() << " ---> "<< energy <<endl;
811 } else {
812 TofData* eTofData = (*it);
813 energy+=eTofData->energy();
814
815 //cout << "Add energy" << endl;
816 }
817 break;
818 }
819 }
820
821 vector<Identifier> NeighborVec=getNextNeighbors(*iSeed);
822 vector<Identifier>::iterator iNeigh;
823 for(iNeigh=NeighborVec.begin(); iNeigh!=NeighborVec.end(); iNeigh++)
824 {
825
826 vector<TofData*>::iterator it2;
827 for(it2=tofDataVec.begin(); it2!=tofDataVec.end(); it2++)
828 {
829 if((*it2)->identify()==*iNeigh) {
830 //cout<<"iNeigh="<<*iNeigh<<endl;
831 if((*it)->barrel()) {
832 TofData* bTofData2 = (*it2);
833
834 if(fabs(bTofData2->zpos())>2) continue;
835 if(m_output) {
836 m_seed_dist = seedPos-bTofData2->zpos();
837 m_tuple2->write();
838 }
839 if(fabs(seedPos-bTofData2->zpos())>0.3) continue;
840 zpos+=bTofData2->zpos()*bTofData2->energy();
841 energy+=bTofData2->energy();
842 //cout << "Add energy barrel neig " << TofID::layer(*iNeigh) <<" " << TofID::phi_module(*iNeigh) << " " << bTofData2->energy() << " ---> "<< energy <<endl;
843 } else {
844 TofData* eTofData2 = (*it2);
845 energy+=eTofData2->energy();
846 }
847 break;
848 }
849 }
850 }
851 if(energy>0) zpos/=energy;
852
853 //for charged track, set energy
854 if(neutral==false) {
855 if(fabs(zpos)<1.15&&energy>5.&&energy<1000) {
856 (*iMatch)->setEnergy(energy/1000);
857 }
858 continue;
859 }
860
861 //for neutral track
862 if(fabs(zpos)<1.15&&energy>5.&&energy<1000) { //shower energy cut = 10MeV
863 RecTofTrack* tof = new RecTofTrack;
864 tof->setTofID(*iSeed);
865 tof->setZrHit(zpos);
866 tof->setEnergy(energy/1000); //MeV-->GeV
867 recTofTrackCol->push_back(tof);
868
869 if(m_output) {
870 m_shower_part = barrel_ec;
871 m_shower_layer = layer;
872 m_shower_im = im;
873 m_shower_zpos = zpos;
874 m_shower_energy = energy;
875 m_tuple1->write();
876 }
877 }
878 }//close if(!is_mrpc)
879 else // mrpc_case
880 {
881 //Determine wether the track is a neutral one or not:(Compare with charged tracks)
882
883 //cout << "TofShower::findShower Seed is mrpc" << endl;
884
885 bool neutral =true;
886 int partId_seed = TofID::barrel_ec(*iSeed);
887 int module_seed = TofID::phi_module(*iSeed)/25;
888 int strip_seed = TofID::phi_module(*iSeed)%25;
889
890
891
892 RecTofTrackCol::iterator rec_charged,iMatch;
893 for( rec_charged = recTofTrackCol->begin(); rec_charged != recTofTrackCol->end();rec_charged++)
894 {
895
896 TofHitStatus *status = new TofHitStatus;
897 status->setStatus((*rec_charged)->status());
898 if((status->is_barrel())) continue; //Consider only endcap tracks
899 if((*rec_charged)->tofID()<=125) continue;
900
901 int partId_charged =-1;
902 int module_charged =-1;
903 int strip_charged =-1;
904
905 if(status->is_east())
906 {
907 if((*rec_charged)->tofID()>=600) {partId_charged=3; module_charged = ((*rec_charged)->tofID()-574)/25; strip_charged =((*rec_charged)->tofID()-574)%25; }
908 else {partId_charged=4; module_charged = ((*rec_charged)->tofID()-100)/25; strip_charged =((*rec_charged)->tofID()-100)%25; }
909 }
910 else
911 {
912 if((*rec_charged)->tofID()>=600) {partId_charged=6; module_charged = ((*rec_charged)->tofID()-574)/25; strip_charged =((*rec_charged)->tofID()-574)%25; }
913 else {partId_charged=5;module_charged = ((*rec_charged)->tofID()-100)/25; strip_charged =((*rec_charged)->tofID()-100)%25;}
914 }
915
916 delete status;
917
918 if( (partId_seed==3 || partId_seed==4) && (partId_charged==3 || partId_charged==4 ) )
919 {
920 if( partId_seed==partId_charged && module_charged == module_seed && abs(strip_charged-strip_seed)<3) {iMatch=rec_charged; neutral =false; break;}
921 if( partId_seed==partId_charged && module_charged == module_seed)
922 {
923 if(strip_seed%2==0 && strip_charged%2!=0 && strip_charged<strip_seed && abs(strip_charged-strip_seed)<4) {iMatch=rec_charged; neutral =false; break;}
924 if(strip_seed%2!=0 && strip_charged%2==0 && strip_charged>strip_seed && abs(strip_charged-strip_seed)<4) {iMatch=rec_charged; neutral =false; break;}
925 }
926
927
928 if( partId_seed!=partId_charged && abs(module_charged-module_seed)<2)
929 {
930 if(strip_seed%2==0 && strip_charged%2!=0 && abs(strip_charged-strip_seed)<3 ) {iMatch=rec_charged; neutral =false; break;}
931 if(strip_seed%2==0 && strip_charged%2!=0 && strip_charged<strip_seed && abs(strip_charged-strip_seed)<4) {iMatch=rec_charged; neutral =false; break;}
932
933 if(strip_seed%2!=0 && strip_charged%2==0 && abs(strip_charged-strip_seed)<3 ) {iMatch=rec_charged; neutral =false; break;}
934 if(strip_seed%2!=0 && strip_charged%2==0 && strip_charged>strip_seed && abs(strip_charged-strip_seed)<4) {iMatch=rec_charged; neutral =false; break;}
935 }
936
937 }
938 else if ( (partId_seed==5 || partId_seed==6) && (partId_charged==5 || partId_charged==6 ) )
939 {
940 if( partId_seed==partId_charged && module_charged ==module_seed && abs(strip_charged-strip_seed)<3) {iMatch=rec_charged; neutral =false; break;}
941
942 if( partId_seed==partId_charged && module_charged == module_seed)
943 {
944 if(strip_seed%2==0 && strip_charged%2!=0 && strip_charged<strip_seed && abs(strip_charged-strip_seed)<4) {iMatch=rec_charged; neutral =false; break;}
945 if(strip_seed%2!=0 && strip_charged%2==0 && strip_charged>strip_seed && abs(strip_charged-strip_seed)<4) {iMatch=rec_charged; neutral =false; break;}
946 }
947
948 if( partId_seed!=partId_charged && abs(module_charged-module_seed)<2)
949 {
950 if(strip_seed%2==0 && strip_charged%2!=0 && abs(strip_charged-strip_seed)<3 ) {iMatch=rec_charged; neutral =false; break;}
951 if(strip_seed%2==0 && strip_charged%2!=0 && strip_charged<strip_seed && abs(strip_charged-strip_seed)<4) {iMatch=rec_charged; neutral =false; break;}
952
953 if(strip_seed%2!=0 && strip_charged%2==0 && abs(strip_charged-strip_seed)<3 ) {iMatch=rec_charged; neutral =false; break;}
954 if(strip_seed%2!=0 && strip_charged%2==0 && strip_charged>strip_seed && abs(strip_charged-strip_seed)<4) {iMatch=rec_charged; neutral =false; break;}
955 }
956 }
957
958 }//Close for rec charegd Track
959
960 //cout << "TofShower::findShower Charged/Neutral ? -> OK" << endl;
961
962
963
964 //Now sum up the energy of the seed and its neighbours:
965
966 double energy=0;
967 double timeinfo_neutral=0;
968 double charge_neutral=0;
969 vector<TofData*>::iterator it;
970 for(it=tofDataVec.begin();it!=tofDataVec.end();it++)
971 {
972 if((*it)->barrel()) continue; //We sum up only endcap data
973
974 if((*it)->identify()==*iSeed)
975 {
976 TofData *eTofData = (*it);
977 energy+=eTofData->energy();
978 timeinfo_neutral= eTofData->tdc();
979 charge_neutral= eTofData->adc();// time over threshold in ns.
980
981 //cout<< "TofShower::findShower : Summing up energy " << endl;
982 //cout << " Seed " << (int)(TofID::phi_module(*iSeed)/25)<< " | " << (int)(TofID::phi_module(*iSeed)%25) << " | " << energy <<endl;
983
984 break;
985 }
986 }
987 //Now add the energy from the neighbours
988
989
990 vector<Identifier> NeighborVec=getMRPC_neighbours(*iSeed);
991 vector<Identifier>::iterator iNeigh;
992 for(iNeigh=NeighborVec.begin(); iNeigh!=NeighborVec.end(); iNeigh++)
993 {
994
995 vector<TofData*>::iterator it2;
996 for(it2=tofDataVec.begin();it2!=tofDataVec.end();it2++)
997 {
998 if((*it2)->identify()==(*iNeigh))
999 {
1000 TofData * eTofData2 = (*it2);
1001 energy+=eTofData2->energy();
1002 //cout << " Neigh " << (int)(TofID::phi_module(*iNeigh)/25)<< " | " << (int)(TofID::phi_module(*iNeigh)%25) << " | " << energy <<endl;
1003 break;
1004 }
1005 }//Close for it2
1006
1007 }//close for i neigh.,
1008
1009
1010
1011 if(neutral==false) // set energy for charged Tracks!
1012 {
1013 if(energy<1000) (*iMatch)->setEnergy(energy/1000.); //MeV->GeV
1014 /*
1015 cout<< "TofShower::findShower --> charged track : Final information" << endl;
1016 cout << "partid_seed " << partId_seed<<endl;
1017 cout << "module_seed " << module_seed<<endl;
1018 cout << "strip_seed " << strip_seed << endl;
1019 cout << "neutral ? " << neutral << endl;
1020 cout << "Pulselength " << charge_neutral <<endl;
1021 cout << "Energy SUM [MeV] " << energy << endl;
1022 cout << "_____________________________________" << endl;
1023 */
1024
1025 continue;
1026 }
1027
1028
1029 if(neutral == true)
1030 {
1031
1032
1033 double tof= timeinfo_neutral-t0;
1034 if(tof<0) tof=0;
1035
1036
1037 RecTofTrack* tof_track =new RecTofTrack;
1038 tof_track->setTofID(*iSeed);
1039 tof_track->setTof(tof);//time of flight in ns
1040 tof_track->setPh(charge_neutral);//collected ADC charge
1041 tof_track->setEnergy(energy/1000.); //MeV->GeV ; Energyloss for the seed and its neighborhood for the whole tof module!!!!
1042 recTofTrackCol->push_back(tof_track);
1043
1044 /*
1045 cout << "TofShower::findShower --> neutral track : Final information" << endl;
1046 cout << "partid_seed " << partId_seed<<endl;
1047 cout << "module_seed " << module_seed<<endl;
1048 cout << "strip_seed " << strip_seed << endl;
1049 cout << "neutral ? " << neutral << endl;
1050 cout << "tof " << timeinfo_neutral-t0 << endl;
1051 cout << "pulselength " << charge_neutral << endl;
1052 cout << "Energy SUM [MeV] " << energy << endl;
1053 cout << "_____________________________________" << endl;
1054 */
1055
1056 }
1057
1058
1059 }//close else mrpc_case
1060
1061 }//close for over seeed
1062}//close find shower!
1063
1064void TofShower::readCalibPar()
1065{
1066 string paraPath = getenv("TOFENERGYRECROOT");
1067 paraPath += "/share/peak.dat";
1068 ifstream in;
1069 in.open(paraPath.c_str());
1070 assert(in);
1071 for(int i=0;i<176;i++) {
1072 in>>m_ecalib[i];
1073 }
1074 in.close();
1075
1076 paraPath = getenv("TOFENERGYRECROOT");
1077 paraPath += "/share/calib.dat";
1078 ifstream in1;
1079 in1.open(paraPath.c_str());
1080 assert(in1);
1081 for(int i=0;i<176;i++) {
1082 in1>>m_calib[i][0]>>m_calib[i][1]>>m_calib[i][2]>>m_calib[i][3];
1083 }
1084 in1.close();
1085}
1086
1087double TofShower::ecalib(const int nsci) const
1088{
1089 if(nsci<176) {
1090 return m_ecalib[nsci];
1091 } else {
1092 return 0;
1093 }
1094}
1095
1096void TofShower::setEcalib(const int nsci, const double ecalib)
1097{
1098 if(nsci<176) {
1099 m_ecalib[nsci]=ecalib;
1100 }
1101}
1102
1103double TofShower::calib(const int n, const int m) const
1104{
1105 if(n<176&&m<4) {
1106 return m_calib[n][m];
1107 } else {
1108 return 0;
1109 }
1110}
1111
1112void TofShower::setCalib(const int n, const int m, const double ecalib)
1113{
1114 if(n<176&&m<4) {
1115 m_calib[n][m]=ecalib;
1116 }
1117}
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128//Be carefull: As I did have problems with the linking, the following function does
1129//exist in TofTrack.cxx in TofRec, too.
1130int TofShower::get_mrpc_partid_neighbour(string which,int stripid,int layer)
1131{
1132
1133 int strip = BesTofDigitizerEcV4::Get_stripnumber_from_unique_identifier(stripid);
1134
1135 int returnvalue=0;
1136
1137 if(which=="l")
1138 {
1139 if(strip%2==0) returnvalue = layer;
1140 else
1141 {
1142 if(layer==3) returnvalue = 4;
1143 if(layer==4) returnvalue = 3;
1144 if(layer==5) returnvalue = 6;
1145 if(layer==6) returnvalue = 5;
1146 }
1147 }
1148
1149 if(which=="r")
1150 {
1151 if(strip%2!=0) returnvalue = layer;
1152 else
1153 {
1154 if(layer==3) returnvalue = 4;
1155 if(layer==4) returnvalue = 3;
1156 if(layer==5) returnvalue = 6;
1157 if(layer==6) returnvalue = 5;
1158 }
1159 }
1160
1161 if(which=="u")returnvalue = layer;
1162 if(which=="b")returnvalue = layer;
1163
1164 return returnvalue;
1165
1166}
1167
1168
1169//Be carefull: As I did have problems with the linking, the following function does
1170//exist in TofTrack.cxx in TofRec, too.
1171int TofShower::get_mrpc_stripid_neighbour(string which,int stripid,int layer)
1172{
1173 int module= BesTofDigitizerEcV4::Get_module_mrpc_from_unique_identifier(stripid);
1174 int strip = BesTofDigitizerEcV4::Get_stripnumber_from_unique_identifier(stripid);
1175
1176 int neighbourstrip=0;
1177 int neighbourmodule=0;
1178
1179 int returnvalue=0;
1180
1181 if(strip==0) return 0; //Falls vorher schon strip 1,2 oder 23,24 gibt die Funtion 0 zurück, diese vlt. nochmalige Eingabe wird hier abgefangen und wieder eine 0 zureuckgegebn!
1182
1183 if(which=="l")
1184 {
1185 if(strip%2==0) //gerade Stripnummer
1186 {
1187 neighbourstrip=strip-1; //linke Nachbar ist dann immer -1
1188 if(layer==3 || layer== 6)
1189 {neighbourmodule=module;}//Nachbarmodul ist selbes Modul
1190 else
1191 {neighbourmodule=module;}//Nachbarmodul ist selbes Modul
1192 }
1193 else //ungerade stripnummer
1194 {
1195 neighbourstrip=strip+1;
1196
1197 if(layer==3 || layer== 5)
1198 {
1199 neighbourmodule=module; //In diesem Fall sind die Nachbarschaftsmodule immer gleich
1200 }
1201 else
1202 {
1203 neighbourmodule=module+1;
1204 if(neighbourmodule==19){neighbourmodule=1;}
1205 }
1206
1207 }
1208
1209 returnvalue= BesTofDigitizerEcV4::Produce_unique_identifier(neighbourmodule, neighbourstrip);
1210 } //close if which==l
1211
1212 else if(which=="r")
1213 {
1214 if(strip%2==0) //gerade Stripnummer
1215 {
1216 neighbourstrip=strip-1; //rechte Nachbar ist dann immer -1
1217
1218 if(layer==4 || layer== 6)
1219 {neighbourmodule=module;}//In diesem Fall ist das rechte Nachbarmodul immer die gleiche modulnummer
1220 else
1221 {neighbourmodule=module-1;
1222 if(neighbourmodule==0){neighbourmodule=18;}
1223 }
1224 }
1225 else //ungerade stripnummer
1226 {
1227 neighbourstrip=strip+1;
1228 if(layer==3 || layer== 6)
1229 {neighbourmodule=module;}//Nachbarmodul ist selbes Modul
1230 else
1231 {neighbourmodule=module;} //Nachbarmodul ist selbes Modul
1232
1233 }
1234
1235 returnvalue=BesTofDigitizerEcV4::Produce_unique_identifier(neighbourmodule, neighbourstrip);//geaendert: neighbourmodule=module in dieser und unter Funktionen entfernt..
1236 } //close if which==r
1237
1238 else if(which=="u")
1239 {
1240 if(strip==23 || strip==24) {neighbourmodule=0; neighbourstrip=0;}
1241 else
1242 {
1243 neighbourstrip=strip+2;
1244 neighbourmodule=module;
1245 }
1246
1247 returnvalue= BesTofDigitizerEcV4::Produce_unique_identifier(neighbourmodule, neighbourstrip);
1248 } //close if which==u
1249
1250 else if(which=="b")
1251 {
1252 if(strip==1 || strip==2) {neighbourmodule=0; neighbourstrip=0;}
1253 else
1254 {
1255 neighbourstrip=strip-2;
1256 neighbourmodule=module;
1257 }
1258
1259 returnvalue=BesTofDigitizerEcV4::Produce_unique_identifier(neighbourmodule, neighbourstrip);
1260 } //close if which==b
1261
1262
1263 return returnvalue;
1264}
1265
1266
1267
n
const Int_t n
Definition: DataBase/tau_mode.c:65
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
energy
************Class m_ypar INTEGER m_KeyWgt INTEGER m_nphot INTEGER m_KeyGPS INTEGER m_IsBeamPolarized INTEGER m_EvtGenInterface DOUBLE PRECISION m_Emin DOUBLE PRECISION m_sphot DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_q2 DOUBLE PRECISION m_PolBeam2 DOUBLE PRECISION m_xErrPb *COMMON c_KK2f $ !CMS energy average $ !Spin Polarization vector first beam $ !Spin Polarization vector second beam $ !Beam energy spread[GeV] $ !minimum hadronization energy[GeV] $ !input READ never touch them !$ !debug facility $ !maximum weight $ !inverse alfaQED $ !minimum real photon energy
Definition: KK2f.h:50
RecTofTrackCol
ObjectVector< RecTofTrack > RecTofTrackCol
Definition: RecTofTrack.h:33
tofQElecSvc
ITofQElecSvc * tofQElecSvc
Definition: TofRawDataProvider.cxx:32
tofQCorrSvc
ITofQCorrSvc * tofQCorrSvc
Definition: TofRawDataProvider.cxx:31
tofCaliSvc
ITofCaliSvc * tofCaliSvc
Definition: TofRawDataProvider.cxx:30
BesTofDigitizerEcV4::Produce_unique_identifier
static G4int Produce_unique_identifier(G4int module_mrpc_f, G4int readoutstripnumber_f)
Definition: BesTofDigitizerEcV4.cc:841
BesTofDigitizerEcV4::Get_module_mrpc_from_unique_identifier
static G4int Get_module_mrpc_from_unique_identifier(G4int unique_identifier_f)
Definition: BesTofDigitizerEcV4.cc:854
BesTofDigitizerEcV4::Get_stripnumber_from_unique_identifier
static G4int Get_stripnumber_from_unique_identifier(G4int unique_identifier_f)
Definition: BesTofDigitizerEcV4.cc:866
DstTofTrack::setEnergy
void setEnergy(double energy)
Definition: DstTofTrack.h:132
DstTofTrack::setPh
void setPh(double ph)
Definition: DstTofTrack.h:96
DstTofTrack::setZrHit
void setZrHit(double zrhit)
Definition: DstTofTrack.h:95
DstTofTrack::setTof
void setTof(double tof)
Definition: DstTofTrack.h:97
DstTofTrack::setTofID
void setTofID(int tofID)
Definition: DstTofTrack.h:91
ITofCaliSvc::BPh
virtual const double BPh(double ADC1, double ADC2, double zHit, unsigned int id)=0
ITofCaliSvc::ZTDC
virtual const double ZTDC(double tleft, double tright, unsigned id)=0
TofData::tofId
int tofId()
Definition: TofData.cxx:494
TofData::setEnergy
void setEnergy(double energy)
Definition: TofData.h:212
TofData::adc2
double adc2()
Definition: TofData.cxx:580
TofData::energy
double energy() const
Definition: TofData.h:197
TofData::setZpos
void setZpos(double zpos)
Definition: TofData.h:211
TofData::adc
double adc()
Definition: TofData.cxx:598
TofData::tdc2
double tdc2()
Definition: TofData.cxx:589
TofData::identify
unsigned int identify() const
Definition: TofData.h:131
TofData::zpos
double zpos() const
Definition: TofData.h:196
TofData::tdc1
double tdc1()
Definition: TofData.cxx:571
TofData::adc1
double adc1()
Definition: TofData.cxx:562
TofData::tdc
double tdc()
Definition: TofData.cxx:607
TofHitStatus::is_barrel
bool is_barrel() const
Definition: TofHitStatus.h:26
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:126
TofHitStatus::is_east
bool is_east() const
Definition: TofHitStatus.h:27
TofID::cell_id
static Identifier cell_id(int barrel_ec, int layer, int phi_module, int end)
For a single crystal.
Definition: TofID.cxx:156
TofID::end
static int end(const Identifier &id)
Definition: TofID.cxx:129
TofID::cell_id_mrpc
static Identifier cell_id_mrpc(int partID, int scinNum)
Definition: TofID.cxx:177
TofID::phi_module
static int phi_module(const Identifier &id)
Definition: TofID.cxx:117
TofID::barrel_ec
static int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition: TofID.cxx:95
TofID::is_mymrpc
static bool is_mymrpc(const Identifier &id)
Definition: TofID.cxx:58
TofID::layer
static int layer(const Identifier &id)
Definition: TofID.cxx:109
TofShower::calib
double calib(const int n, const int m) const
Definition: TofShower.cxx:1103
TofShower::ecalib
double ecalib(const int nsci) const
Definition: TofShower.cxx:1087
TofShower::readCalibPar
void readCalibPar()
Definition: TofShower.cxx:1064
TofShower::findSeed
void findSeed(vector< TofData * > &tofDataVec)
Definition: TofShower.cxx:626
TofShower::BookNtuple
void BookNtuple(NTuple::Tuple *&tuple, NTuple::Tuple *&tuple1, NTuple::Tuple *&tuple2)
Definition: TofShower.cxx:33
TofShower::setEcalib
void setEcalib(const int nsci, const double ecalib)
Definition: TofShower.cxx:1096
TofShower::get_mrpc_stripid_neighbour
int get_mrpc_stripid_neighbour(string which, int stripid, int layer)
Definition: TofShower.cxx:1171
TofShower::getMRPC_neighbours
vector< Identifier > getMRPC_neighbours(const Identifier &id)
Definition: TofShower.cxx:528
TofShower::energyCalib
void energyCalib(vector< TofData * > &tofDataVec, RecTofTrackCol *recTofTrackCol)
Definition: TofShower.cxx:71
TofShower::getNextNeighbors
vector< Identifier > getNextNeighbors(const Identifier &id)
Definition: TofShower.cxx:471
TofShower::setCalib
void setCalib(const int n, const int m, const double ecalib)
Definition: TofShower.cxx:1112
TofShower::getNeighbors
vector< Identifier > getNeighbors(const Identifier &id)
Definition: TofShower.cxx:391
TofShower::findShower
void findShower(vector< TofData * > &tofDataVec, RecTofTrackCol *recTofTrackCol, double)
Definition: TofShower.cxx:745
TofShower::get_mrpc_partid_neighbour
int get_mrpc_partid_neighbour(string which, int stripid, int layer)
Definition: TofShower.cxx:1130
num
int num[96]
Definition: ranlxd.c:373