10 #include <TLorentzVector.h>
28 #include "TGraphErrors.h"
100 vector< vector< double > >
z_track_graph(768, vector< double >() );
101 vector< vector< double > >
z_point_graph(768, vector< double >() );
102 vector< vector< double > >
z_deltat_graph(768, vector< double >() );
104 vector< vector< double > >
z_track_thrown(768, vector< double >() );
105 vector< vector< double > >
z_point_thrown(768, vector< double >() );
108 #include <JANA/JApplication.h>
109 #include <JANA/JFactory.h>
146 gPARMS->SetDefaultParameter(
"BCALPOINTCALIB:DEBUG",
DEBUG,
"Control the creation of extra histograms");
147 gPARMS->SetDefaultParameter(
"BCALPOINTCALIB:DEBUG",
VERBOSE,
"Control verbose output");
154 japp->RootWriteLock();
159 TDirectory *
main = gDirectory;
160 gDirectory->mkdir(
"bcal_point_calibs")->cd();
167 h1_Num_matched_showers =
new TH1I(
"h1_Num_matched_showers",
"Number of matched showers per event", 20,0,20);
176 h1_thrown_per_channel =
new TH1D(
"h1_thrown_per_channel",
"Percentage of thrown points per Channel",800,0,800);
181 h1pt_Num_points =
new TH1I(
"h1pt_Num_points",
"pt points per shower",20,0,20);
188 h1_E =
new TH1I(
"h1_E",
"Matched energy per shower",nbins,0,2);
189 h1_E->SetXTitle(
"Energy per shower (GeV)");
190 h1_E->SetYTitle(
"counts");
191 h1_E_raw =
new TH1I(
"h1_E_raw",
"Matched raw energy per shower",nbins,0,2);
192 h1_E_raw->SetXTitle(
"Raw Energy per shower (GeV)");
194 h1_x =
new TH1I(
"h1_x",
"x pos of shower",nbins2,-100,100);
195 h1_x->SetXTitle(
"x (cm)");
196 h1_x->SetYTitle(
"counts");
197 h1_y =
new TH1I(
"h1_y",
"y pos of shower",nbins2,-100,100);
198 h1_y->SetXTitle(
"y (cm)");
199 h1_y->SetYTitle(
"counts");
200 h1_z =
new TH1I(
"h1_z",
"z pos of shower",nbins2,0,500);
201 h1_z->SetXTitle(
"z (cm)");
202 h1_z->SetYTitle(
"counts");
203 h1_t =
new TH1I(
"h1_t",
"time of shower",nbins2,-100,100);
204 h1_t->SetXTitle(
"time of shower (ns)");
205 h1_t->SetYTitle(
"counts");
206 h1_N_cell =
new TH1I(
"h1_N_cell",
"N_cell of shower",20,0,20);
211 h1trk_FOM =
new TH1I(
"h1trk_FOM",
"FOM for matched tracks", nbins,0,1);
214 h1trk_pmag =
new TH1I(
"h1trk_mag",
"pmag for matched tracks", nbins,0,4);
217 h1trk_px =
new TH1I(
"h1trk_px",
"px for matched tracks", nbins,0,2);
220 h1trk_py =
new TH1I(
"h1trk_py",
"py for matched tracks", nbins,0,2);
221 h1trk_py->SetXTitle(
"Number of tracks");
223 h1trk_pz =
new TH1I(
"h1trk_pz",
"pz for matched tracks", nbins,0,4);
226 h1trk_x =
new TH1I(
"h1trk_x",
"x for matched tracks", nbins2,-30,30);
229 h1trk_y =
new TH1I(
"h1trk_y",
"y for matched tracks", nbins2,-30,30);
232 h1trk_z =
new TH1I(
"h1trk_z",
"z for matched tracks", nbins2,-50,250);
235 h1trk_energy =
new TH1I(
"h1trk_energy",
"energy for matched tracks", nbins,0,4);
239 h2_Evsenergy =
new TH2I(
"h2_Evsenergy",
"E vs energy matched tracks", nbins,0,4,nbins,0,4);
242 h2_pmagvsenergy =
new TH2I(
"h2_pmagvsenergy",
"pmag vs energy matched tracks", nbins,0,4,nbins,0,4);
247 h1pt_module =
new TH1I(
"h1pt_module",
"pt module number",50,0,50);
250 h1pt_layer =
new TH1I(
"h1pt_layer",
"pt layer number",5,0,5);
253 h1pt_sector =
new TH1I(
"h1pt_sector",
"pt sector number",5,0,5);
256 h1pt_cell_id =
new TH1I(
"h1pt_cell_id",
"pt cell id number",800,0,800);
259 h1pt_energy =
new TH1I(
"h1pt_energy",
"pt energy",nbins,0,1);
262 h1pt_energy_US =
new TH1I(
"h1pt_energy_US",
"pt energy US",nbins,0,1);
265 h1pt_energy_DS =
new TH1I(
"h1pt_energy_DS",
"pt energy DS",nbins,0,1);
268 h1pt_z =
new TH1I(
"h1pt_z",
"z for point",nbins2,0,500);
269 h1pt_z->SetXTitle(
"z (cm)");
270 h1pt_z->SetYTitle(
"counts");
271 h1pt_r =
new TH1I(
"h1pt_r",
"r for point",nbins,60,90);
272 h1pt_r->SetXTitle(
"r (cm)");
273 h1pt_r->SetYTitle(
"counts");
274 h1pt_r_size =
new TH1I(
"h1pt_r_size",
"r size of point", 20, 0, 20);
277 h1pt_phi =
new TH1I(
"h1pt_phi",
"phi of point", 700, 0, 7);
280 h1pt_t =
new TH1I(
"h1pt_t",
"t for point",nbins2,-100,100);
282 h1pt_t->SetYTitle(
"counts");
295 h2_zpoint_vs_ztrack =
new TH2I(
"h2_zpoint_vs_ztrack",
"z point vs z track for the BCAL",1000,0,500,1000,0,500);
298 h2_zpoint_vs_ztrack_thrown =
new TH2I(
"h2_zpoint_vs_ztrack_thrown",
"thrown z point vs z track for the BCAL",1000,0,500,1000,0,500);
322 Run_Number = runnumber;
351 loop->GetSingle(locTrigger);
356 vector<const DBCALGeometry *> BCALGeomVec;
357 loop->Get(BCALGeomVec);
358 if(BCALGeomVec.size() == 0)
359 throw JException(
"Could not load DBCALGeometry object!");
360 auto locBCALGeom = BCALGeomVec[0];
362 vector<const DTrackFitter *> fitters;
365 if(fitters.size()<1){
366 _DBG_<<
"Unable to get a DTrackFinder object!"<<endl;
367 return RESOURCE_UNAVAILABLE;
372 vector<const DBCALShower*> locBCALShowers;
373 vector<const DBCALHit*> bcalhits;
374 vector<const DBCALPoint*> locBCALPoints;
376 loop->Get(locBCALShowers);
378 loop->Get(locBCALPoints);
380 vector<const DTrackTimeBased*> locTrackTimeBased;
381 loop->Get(locTrackTimeBased);
383 vector <const DBCALShower *> matchedShowers;
384 vector < const DBCALShower *> matchedShowersneg;
385 vector < const DBCALShower *> matchedShowerspos;
386 vector <const DTrackTimeBased* > matchedTracks;
389 map< const DBCALShower*, vector<const DBCALPoint*> > matchedShowerPoints_cache;
390 map< const DBCALPoint *, const DBCALUnifiedHit *> upstreamHit_cache;
391 map< const DBCALPoint *, const DBCALUnifiedHit *> downstreamHit_cache;
395 for (
unsigned int i=0; i < locTrackTimeBased.size() ; ++i){
396 vector<DTrackFitter::Extrapolation_t>extrapolations=locTrackTimeBased[i]->extrapolations.at(
SYS_BCAL);
397 if (extrapolations.size()==0)
continue;
399 for (
unsigned int j=0; j< locBCALShowers.size(); ++j){
401 double x = locBCALShowers[j]->x;
402 double y = locBCALShowers[j]->y;
403 double z = locBCALShowers[j]->z;
406 double R = pos_bcal.Perp();
407 double phi = pos_bcal.Phi();
410 double dPhi = mypos.Phi()-phi;
411 if (dPhi<-M_PI) dPhi+=2.*M_PI;
412 if (dPhi>M_PI) dPhi-=2.*M_PI;
414 double dZ = TMath::Abs(mypos.Z() - z);
417 if(dZ < 30.0 && fabs(dPhi) < 0.18) {
418 matchedShowers.push_back(locBCALShowers[j]);
419 matchedTracks.push_back(locTrackTimeBased[i]);
437 const float*bcal_radii = locBCALGeom->GetBCAL_radii();
440 double R1 = 0.5*(bcal_radii[0] + bcal_radii[1]);
441 double R2 = 0.5*(bcal_radii[1] + bcal_radii[2]);
442 double R3 = 0.5*(bcal_radii[2] + bcal_radii[3]);
443 double R4 = 0.5*(bcal_radii[3] + bcal_radii[4]);
446 Int_t numshowers_per_event = matchedShowers.size();
447 if (numshowers_per_event > 0)
449 Int_t numtracks_per_event = matchedTracks.size();
450 if (numtracks_per_event > 0)
454 for(
int i=0; i<numshowers_per_event; i++) {
455 vector<const DBCALPoint*> points;
456 matchedShowers[i]->Get(points);
457 matchedShowerPoints_cache[ matchedShowers[i] ] = points;
459 vector<const DBCALUnifiedHit*> hits;
460 for(
unsigned int k=0; k<points.size(); k++) {
461 points[k]->Get(hits);
467 upstreamHit_cache[ points[k] ] = hits[0];
468 downstreamHit_cache[ points[k] ] = hits[1];
470 upstreamHit_cache[ points[k] ] = hits[1];
471 downstreamHit_cache[ points[k] ] = hits[0];
478 japp->RootFillLock(
this);
480 for(
int i=0; i<numshowers_per_event; i++){
483 Float_t E = matchedShowers[i]->E;
484 Float_t E_raw = matchedShowers[i]->E_raw;
485 Float_t
x = matchedShowers[i]->x;
486 Float_t
y = matchedShowers[i]->y;
487 Float_t z = matchedShowers[i]->z;
488 Float_t t = matchedShowers[i]->t;
489 Int_t N_cell = matchedShowers[i]->N_cell;
500 double FOM = matchedTracks[i]->FOM;
501 double pmag = matchedTracks[i]->pmag();
502 double px = matchedTracks[i]->px();
503 double py = matchedTracks[i]->py();
504 double pz = matchedTracks[i]->pz();
505 double x_track = matchedTracks[i]->x();
506 double y_track = matchedTracks[i]->y();
507 double z_track = matchedTracks[i]->z();
509 double energy = matchedTracks[i]->energy();
525 vector<DTrackFitter::Extrapolation_t>extrapolations=matchedTracks[i]->extrapolations.at(
SYS_BCAL);
532 double mypos_z_array[4] = {};
533 mypos_z_array[0] = mypos_1.Z();
534 mypos_z_array[1] = mypos_2.Z();
535 mypos_z_array[2] = mypos_3.Z();
536 mypos_z_array[3] = mypos_4.Z();
538 double target_center = 65.0;
544 vector<const DBCALPoint*> &points = matchedShowerPoints_cache[ matchedShowers[i] ];
547 Int_t numpoints_per_shower = points.size();
550 for (
int jj=0; jj<numpoints_per_shower; jj++) {
551 int module = points[jj]->module();
552 int layer = points[jj]->layer();
553 int sector = points[jj]->sector();
554 int cell_id = (module-1)*16 + (layer-1)*4 + sector;
555 double energy = points[jj]->E();
556 double energy_US = points[jj]->E_US();
557 double energy_DS = points[jj]->E_DS();
558 double z_point_wrt_target_center = points[jj]->z();
559 double z_point = z_point_wrt_target_center + target_center;
560 double r_point = points[jj]->r();
561 double phi_point = points[jj]->phi();
562 double t_point = points[jj]->t();
567 double deltat_point = up_hit->
t_ADC - down_hit->
t_ADC;
586 if(mypos_z_array[layer-1] > 0 && z_point > 0){
597 if(TMath::Abs(mypos_z_array[layer-1] - z_point) < 400.0){
616 japp->RootFillUnLock(
this);
644 japp->RootFillLock(
this);
649 TDirectory *
main = gDirectory;
650 TDirectory *
locDirectory = (TDirectory*)gDirectory->FindObjectAny(
"bcal_point_calibs");
658 double percent[768] = {0};
659 double number_of_points = 0.0;
660 double number_of_thrown_points = 0.0;
663 for(
int m = 0; m < channels; ++m){
665 h2_tgraph[m]->SetTitle(Form(
"z of point vs z of track for channel %i", m+1));
666 h2_tgraph[m]->GetXaxis()->SetTitle(
"z of track (cm)");
667 h2_tgraph[m]->GetYaxis()->SetTitle(
"z of point (cm)");
668 h2_tgraph[m]->GetXaxis()->SetLimits(0.0, 500.0);
669 h2_tgraph[m]->GetYaxis()->SetRangeUser(0.0, 500.0);
673 h2_tgraph[m]->Write(Form(
"h2_tgraph[%i]", m+1));
676 h2_tgraph_deltat[m]->SetTitle(Form(
"Delta T of point vs z of track for channel %i", m+1));
688 h2_tgraph_thrown[m]->SetTitle(Form(
"Thrown z of point vs z of track for channel %i", m+1));
700 percent[m] = 100 * (number_of_thrown_points / number_of_points);
701 if(percent[m] >= 0.0 && percent[m] <= 100.0){
708 japp->RootFillUnLock(
this);
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
bool ExtrapolateToRadius(double R, const vector< Extrapolation_t > &extraps, DVector3 &pos, DVector3 &mom, double &t, double &s) const
vector< vector< double > > z_track_thrown(768, vector< double >())
static TH1I * h1pt_sector
static TH1I * h1_ztrack_minus_zpoint
static TH1D * h1_thrown_per_channel
static TH1I * h1pt_r_size
~JEventProcessor_BCAL_point_calib()
vector< vector< double > > z_track_graph(768, vector< double >())
static TH1I * h1pt_energy_DS
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
uint32_t Get_L1FrontPanelTriggerBits(void) const
vector< vector< double > > z_deltat_graph(768, vector< double >())
static TH1I * h1pt_module
static TGraph * h2_tgraph_thrown[768]
static TH2I * h2_pmagvsenergy
static TH2I * h2_zpoint_vs_ztrack_thrown
TDirectory * locDirectory
static TH1I * h1pt_energy
static TH2I * h2_Evsenergy
JEventProcessor_BCAL_point_calib()
static TH1I * h1pt_Num_points
vector< vector< double > > z_point_graph(768, vector< double >())
jerror_t init(void)
Called once at program start.
float t_ADC
Time from fADC.
static TGraph * h2_tgraph[768]
static TH1I * h1pt_energy_US
static TH1I * h1trk_Num_matched_tracks
static TH1I * h1trk_energy
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
static TH2I * h2_zpoint_vs_ztrack
const DTrackFitter * fitter
jerror_t fini(void)
Called after last event of last event source has been processed.
static TGraph * h2_tgraph_deltat[768]
printf("string=%s", string)
static TH1I * h1_abs_ztrack_minus_zpoint
vector< vector< double > > z_point_thrown(768, vector< double >())
static TH2I * h2_ztrack_minus_zpoint_vs_ztrack
int main(int argc, char *argv[])
static TH1I * h1_Num_matched_showers
static TH1I * h1pt_cell_id