10 #include <TLorentzVector.h>
35 InitJANAPlugin(locApplication);
40 #define FCAL_Z_OFFSET 640.0-65.0
46 TDirectory *
dir =
new TDirectoryFile(
"BCAL",
"BCAL");
49 two_gamma_mass =
new TH1F(
"two_gamma_mass",
"two_gamma_mass",2000, 0.0, 3.00);
58 mass_summed =
new TH1F(
"inv_mass",
"inv_mass",100,0.0,1.0);
60 z_shower =
new TH1F(
"z_shower",
"z_shower",1600, 0.0, 6.0);
61 E_shower =
new TH1F(
"E_shower",
"E_shower", 1600, 0.0, 6.0);
62 Neutral_E =
new TH1F(
"Neutral_E",
"Neutral_E",1000,0.0, 400.0);
63 VertexZ =
new TH1F(
"VertexZ",
"VertexZ",1000,-200,300);
64 goodVertexZ =
new TH1F(
"goodVertexZ",
"goodVertexZ",1000,-200,300);
74 E1_Vs_E2CosTheta =
new TH2F(
"E1 Vs E2(1-cosTheta)",
"E1 Vs E2(1-cosTheta)",2000,0.0,1.5,2000.0,0.0,5.0);
75 Thrown_Gamma_Theta =
new TH1F(
"Thrown Gamma Theta",
"Thrown Gamma Theta",300,0.0,150.0);
76 Thrown_Inv_Mass =
new TH1F(
"Thrown Gamma Inv Mass",
"Thrown Gamma Inv Mass",1000,0.0,1.50);
77 Thrown_Vertex =
new TH1F(
"Thrown Gamma Vertex",
"Thrown Gamma Vertex",2000,-300.0,400.0);
81 Point_Energy =
new TH1F(
"Point Energy",
"Point Energy",1000,0.0,2.0);
83 Point_z =
new TH1F(
"Point z",
"Point z",1000,-100.0,600.0);
84 Assoc_E =
new TH1F(
"Assoc E",
"Assoc E",1000,0,2.0);
85 Layer1_Energy_vs_Channel =
new TH2F(
"Layer1 Energy vs channel",
"Layer1 Energy vs channel",768,0.0,768.0,1000,0.0,0.5);
86 Layer2_Energy_vs_Channel =
new TH2F(
"Layer2 Energy vs channel",
"Layer2 Energy vs channel",768,0.0,768.0,1000,0.0,0.5);
87 Layer3_Energy_vs_Channel =
new TH2F(
"Layer3 Energy vs channel",
"Layer3 Energy vs channel",768,0.0,768.0,1000,0.0,0.5);
88 Layer4_Energy_vs_Channel =
new TH2F(
"Layer4 Energy vs channel",
"Layer4 Energy vs channel",768,0.0,768.0,1000,0.0,0.5);
89 Layer1_Energy =
new TH1F(
"Layer1 Energy",
"Layer1 Energy",1500,0.0,1.5);
90 Layer2_Energy =
new TH1F(
"Layer2 Energy",
"Layer2 Energy",1500,0.0,1.5);
91 Layer3_Energy =
new TH1F(
"Layer3 Energy",
"Layer3 Energy",1500,0.0,1.5);
92 Layer4_Energy =
new TH1F(
"Layer4 Energy",
"Layer4 Energy",1500,0.0,1.5);
93 Layer1_Energy_v2 =
new TH1F(
"Layer1 Energy v2",
"Layer1 Energy v2",1500,0.0,1.5);
94 Layer2_Energy_v2 =
new TH1F(
"Layer2 Energy v2",
"Layer2 Energy v2",1500,0.0,1.5);
95 Layer3_Energy_v2 =
new TH1F(
"Layer3 Energy v2",
"Layer3 Energy v2",1500,0.0,1.5);
96 Layer4_Energy_v2 =
new TH1F(
"Layer4 Energy v2",
"Layer4 Energy v2",1500,0.0,1.5);
97 Layer1_Energy_pos =
new TH1F(
"Layer1 Energy pos",
"Layer1 Energy pos",1500,0.0,0.1);
98 Layer2_Energy_pos =
new TH1F(
"Layer2 Energy pos",
"Layer2 Energy pos",1500,0.0,0.1);
99 Layer3_Energy_pos =
new TH1F(
"Layer3 Energy pos",
"Layer3 Energy pos",1500,0.0,0.15);
100 Layer4_Energy_pos =
new TH1F(
"Layer4 Energy pos",
"Layer4 Energy pos",1500,0.0,0.15);
105 Time_Diff =
new TH1F(
"Time Diff",
"Time Diff",3000,0.0,500);
106 Eoverp_plus_nocuts =
new TH1F(
" E over p plus no cuts ",
" E over p plus no cuts ", 1000,0.0,5.0);
107 Eoverp_minus_nocuts =
new TH1F(
" E over p minus no cuts ",
" E over p minus no cuts ",1000,0.0,5.0);
108 Eoverp_plus_cuts =
new TH1F(
" E over p plus cuts ",
" E over p plus cuts ", 1000,0.0,5.0);
109 Eoverp_minus_cuts =
new TH1F(
" E over p minus cuts ",
" E over p minus cuts ",1000,0.0,5.0);
110 Evsp_plus =
new TH2F(
" E vs p plus ",
" E vs p plus", 1000,0.0,10.0, 1000, 0.0, 10.0);
111 Evsp_minus =
new TH2F(
" E vs p minus ",
" E vs p minus", 1000,0.0,10.0, 1000, 0.0, 10.0);
112 Eoverpvsp_plus =
new TH2F(
" E over p vs p plus " ,
" E over p vs p plus " , 1000, 0.0, 10.0, 1000, 0.0, 10.0);
113 Eoverpvsp_minus =
new TH2F(
" E over p vs p minus " ,
" E over p vs p minus " , 1000,0.0,10.0,1000,0.0,10.0);
114 FCAL_Evsp =
new TH2F(
" FCAL E vs p " ,
" FCAL E vs p " ,500,0.0,5.0,250,0.0,2.5);
115 FCAL_Eoverpvsp =
new TH2F(
" FCAL Eoverp vs p " ,
" FCAL Eoverp vs p " , 500, 0.0, 5.0, 120, 0.0, 1.2);
116 FCAL_Eoverp_cuts =
new TH1F (
" FCAL E over p cuts " ,
" FCAL E over p cuts " , 1000, 0.0, 1.0);
117 FCAL_Eoverp_nocuts =
new TH1F(
" FCAL E over p no cuts " ,
" FCAL E over p no cuts " , 1000,0.0,1.0);
120 All_Layer_Energy =
new TH1F(
"All Layer Point Energy",
"All Layer Point Energy",1500,0.0,1.0);
127 Etot_hits =
new TH1F(
"Etot_hits",
"Sum of all hits (GeV)", 200, 0.0, 6.0);
402 vector<const DMCThrown*> locMCThrowns;
403 locEventLoop->Get(locMCThrowns);
405 vector<const DBCALShower*> locBCALShowers;
406 vector<const DFCALShower*> locFCALShowers;
407 vector<const DBCALTruthShower*> truthshowers;
408 vector<const DMCThrown*> mcthrowns;
409 vector<const DBCALHit*> bcalhits;
410 vector<const DBCALCluster*> locBCALClusters;
411 vector<const DFCALCluster*> locFCALClusters;
412 vector<const DBCALPoint*> locBCALPoints;
413 vector<const DVertex*> kinfitVertex;
415 locEventLoop->GetSingle(locDetectorMatches);
416 locEventLoop->Get(locBCALShowers);
417 locEventLoop->Get(locFCALShowers);
418 locEventLoop->Get(bcalhits);
419 locEventLoop->Get(mcthrowns);
420 locEventLoop->Get(truthshowers);
421 locEventLoop->Get(locBCALClusters);
422 locEventLoop->Get(locFCALClusters);
423 locEventLoop->Get(locBCALPoints);
424 locEventLoop->Get(kinfitVertex);
426 vector<const DTrackTimeBased*> locTrackTimeBased;
427 locEventLoop->Get(locTrackTimeBased);
434 vector <const DChargedTrackHypothesis*> locChargedTrackHypothesis;
435 locEventLoop->Get(locChargedTrackHypothesis);
452 vector <const DChargedTrackHypothesis*> locParticles;
453 locEventLoop->Get(locParticles);
458 japp->RootWriteLock();
460 vector <const DBCALShower *> matchedShowers;
461 vector < const DBCALShower *> matchedShowersneg;
462 vector < const DBCALShower *> matchedShowerspos;
463 vector <const DTrackTimeBased* > matchedTracks;
464 vector <const DFCALShower *> matchedFCALShowers;
465 vector < const DFCALCluster *> matchedFCALClusters;
476 for (
unsigned int i=0; i < locTrackTimeBased.size() ; ++i){
477 for (
unsigned int j=0; j< locBCALShowers.size(); ++j){
479 double x = locBCALShowers[j]->x;
480 double y = locBCALShowers[j]->y;
481 double z = locBCALShowers[j]->z;
482 double E = locBCALShowers[j]->E;
484 double R = pos_bcal.Perp();
485 double phi = pos_bcal.Phi();
486 double L2 = 0.81*2.54+65.0;
487 double L3 = L2 + 0.81*2.54*2;
488 double L4 = L3 + 0.81*2.54*3;
489 double L5 = L4 + 0.97*2.54*4;
490 double L41 = L4 + 0.97*2.54*4*1/4;
491 double L42 = L4 + 0.97*2.54*4*2/4;
492 double L43 = L4 + 0.97*2.54*4*3/4;
493 double FOM = TMath::Prob(locTrackTimeBased[i]->chisq, locTrackTimeBased[i]->Ndof);
494 locTrackTimeBased[i]->rt->GetIntersectionWithRadius(R, mypos);
495 locTrackTimeBased[i]->rt->GetIntersectionWithRadius(65.0,myposL1);
496 locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L2,myposL2);
497 locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L3,myposL3);
498 locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L4,myposL4);
499 locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L5,myposL5);
500 locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L41,myposL41);
501 locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L42,myposL42);
502 locTrackTimeBased[i]->rt->GetIntersectionWithRadius(L43,myposL43);
504 double q = locTrackTimeBased[i]->rt->q;
505 p = locTrackTimeBased[i]->momentum().Mag();
506 double dPhi = TMath::Abs(mypos.Phi()-pos_bcal.Phi());
507 double dZ = TMath::Abs(mypos.Z() -
z);
509 if( mypos.Perp() == R && p > 0.8 && dZ < 30.0 && dPhi < 0.18 && q == 1.0)
Eoverp_plus_cuts->Fill(E/p);
510 if( mypos.Perp() == R && p > 0.8 && dZ < 30.0 && dPhi < 0.18 && q == -1.0)
Eoverp_minus_cuts->Fill(E/p);
511 if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==1.0)
Eoverp_plus_nocuts->Fill(E/p);
512 if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==-1.0)
Eoverp_minus_nocuts->Fill(E/p);
513 if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==1.0)
Evsp_plus->Fill(p,E);
514 if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==-1.0)
Evsp_minus->Fill(p,E);
515 if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==1.0)
Eoverpvsp_plus->Fill(p,E/p);
516 if( mypos.Perp() == R && dZ < 30.0 && dPhi < 0.18 && q ==-1.0)
Eoverpvsp_minus->Fill(p,E/p);
518 if(dZ < 30.0 && dPhi < 0.18 && mypos.Perp() == R) {
519 matchedShowers.push_back(locBCALShowers[j]);
520 matchedTracks.push_back(locTrackTimeBased[i]);
524 if(dZ < 30.0 && dPhi < 0.18 && mypos.Perp() == R && q == -1.0){
525 matchedShowersneg.push_back(locBCALShowers[j]);
527 if(dZ < 30.0 && dPhi < 0.18 && mypos.Perp() == R && q == 1.0){
528 matchedShowerspos.push_back(locBCALShowers[j]);
533 for(
unsigned int i = 0 ; i < locTrackTimeBased.size() ; ++i)
535 for(
unsigned int j = 0 ; j < locFCALShowers.size(); ++j)
546 locTrackTimeBased[i]->rt->GetIntersectionWithPlane(fcalpos,norm,pos);
548 double diffX = TMath::Abs(x - pos.X());
549 double diffY = TMath::Abs(y - pos.Y());
550 if(diffX < 3.0 && diffY < 3.0)
552 matchedFCALShowers.push_back(locFCALShowers[j]);
558 for(
unsigned int i = 0 ; i < locTrackTimeBased.size(); ++i)
560 for(
unsigned int j = 0 ; j < locFCALClusters.size(); ++j)
570 locTrackTimeBased[i]->rt->GetIntersectionWithPlane(fcalpos,norm,pos);
572 double diffX = TMath::Abs(x - pos.X());
573 double diffY = TMath::Abs(y - pos.Y());
574 if(diffX < 3.0 && diffY < 3.0)
576 matchedFCALClusters.push_back(locFCALClusters[j]);
590 double locVertexZ, locVertexY, locVertexX;
594 if(locParticles.size() == 0)
597 double locDOCA, locDOCA2, locSmallestDOCA, locTime0;
598 double locAverageTime = 0.0;
604 locSmallestDOCA = 9.9E9;
605 if(locParticles.size()>1){
607 for(
int j = 0; j < (int(locParticles.size())-1); ++j)
610 for(
size_t k= j+1; k < locParticles.size(); ++k)
614 locSmallestDOCA = locDOCA;
615 locVertex = locTempVertex;
616 locVertexZ = locVertex.Z();
617 locVertexY = locVertex.Y();
618 locVertexX = locVertex.X();
624 double kinfitVertexX, kinfitVertexY, kinfitVertexZ, kinfitVertexT;
625 for (
int i = 0 ; i < kinfitVertex.size(); i++)
627 kinfitVertexX = kinfitVertex[i]->dSpacetimeVertex.X();
628 kinfitVertexY = kinfitVertex[i]->dSpacetimeVertex.Y();
629 kinfitVertexZ = kinfitVertex[i]->dSpacetimeVertex.Z();
630 kinfitVertexT = kinfitVertex[i]->dSpacetimeVertex.T();
658 double Etot_reconstructed = 0.0;
659 for(
unsigned int i=0; i<locBCALShowers.size(); i++)
668 Etot_reconstructed += s->
E;
675 for(
unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
679 if (find(matchedShowers.begin(), matchedShowers.end(),locBCALShowers[loc_i]) != matchedShowers.end())
continue;
681 const DBCALShower* a_shower = locBCALShowers[loc_i];
682 vector<const DBCALCluster*> associated_clusters;
683 a_shower->Get(associated_clusters);
685 double z_one = s1->
z;
689 double E_shower_raw = a_shower->
E_raw;
693 for(
unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
696 const DBCALCluster* a_cluster = associated_clusters[loc_j];
697 vector<const DBCALPoint*> associated_points;
698 a_cluster->Get(associated_points);
700 if(associated_points.size() == 0) cout <<
" assoc points is empty " << endl;
702 for(
unsigned int loc_k = 0; loc_k < associated_points.size(); loc_k++)
704 int module2 = associated_points[loc_k]->module();
705 int layer = associated_points[loc_k]->layer();
706 double point_energy = associated_points[loc_k]->E();
726 for(
unsigned int i=0; i< locBCALPoints.size(); i++)
742 for(
unsigned int i=0; i< locBCALClusters.size(); i++)
852 for(
unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
854 const DBCALShower* locBCALShower = locBCALShowers[loc_i];
855 double E_before = locBCALShower->
E;
857 if(find(matchedShowers.begin(), matchedShowers.end(), locBCALShower) == matchedShowers.end())
continue;
859 double x = locBCALShower->
x;
860 double y = locBCALShower->
y;
861 double z = locBCALShower->
z;
862 double E1 = locBCALShower->
E;
863 double first_pos = TMath::Sqrt( TMath::Power(myposL1.X(),2)+TMath::Power(myposL1.Y(),2)+TMath::Power(myposL1.Z(),2));
864 double second_pos = TMath::Sqrt( TMath::Power(myposL2.X(),2)+TMath::Power(myposL2.Y(),2)+TMath::Power(myposL2.Z(),2));
865 double third_pos = TMath::Sqrt( TMath::Power(myposL3.X(),2)+TMath::Power(myposL3.Y(),2)+TMath::Power(myposL3.Z(),2));
866 double fourth_pos = TMath::Sqrt( TMath::Power(myposL4.X(),2)+TMath::Power(myposL4.Y(),2)+TMath::Power(myposL4.Z(),2));
867 double fifth_pos = TMath::Sqrt( TMath::Power(myposL5.X(),2)+TMath::Power(myposL5.Y(),2)+TMath::Power(myposL5.Z(),2));
868 double fourth_pos2_test = TMath::Sqrt( TMath::Power(myposL4.X()-myposL41.X(),2)+TMath::Power(myposL4.Y()-myposL41.Y(),2)+TMath::Power(myposL4.Z()-myposL41.Z(),2));
869 double fourth_pos3_test = TMath::Sqrt( TMath::Power(myposL41.X()-myposL42.X(),2)+TMath::Power(myposL41.Y()-myposL42.Y(),2)+TMath::Power(myposL41.Z()-myposL42.Z(),2));
870 double fourth_pos4_test = TMath::Sqrt( TMath::Power(myposL42.X()-myposL43.X(),2)+TMath::Power(myposL42.Y()-myposL43.Y(),2)+TMath::Power(myposL42.Z()-myposL43.Z(),2));
871 double fourth_pos5_test = TMath::Sqrt( TMath::Power(myposL43.X()-myposL5.X(),2)+TMath::Power(myposL43.Y()-myposL5.Y(),2)+TMath::Power(myposL43.Z()-myposL5.Z(),2));
873 double less_approx_dist = fourth_pos2_test + fourth_pos3_test + fourth_pos4_test + fourth_pos5_test;
884 double E_shower_raw = locBCALShower->
E_raw;
885 double x1 = locBCALShower->
x - kinfitVertexX;
886 double y1 = locBCALShower->
y - kinfitVertexY;
887 double z1 = locBCALShower->
z - kinfitVertexZ;
888 double t1 = locBCALShower->
t;
889 double R1 =
sqrt(x1*x1 + y1*y1 + z1*z1);
891 TLorentzVector
p1(E1*x1/R1, E1*y1/R1, E1*z1/R1, E1);
892 vector<const DBCALCluster*> associated_clusters;
893 locBCALShower->Get(associated_clusters);
895 if(associated_clusters.size() == 0) cout <<
" assoc cluster is empty " << endl;
897 for(
unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
900 const DBCALCluster* a_cluster = associated_clusters[loc_j];
901 vector<const DBCALPoint*> associated_points;
902 a_cluster->Get(associated_points);
904 if(associated_points.size() == 0) cout <<
" assoc points is empty " << endl;
905 double Layer1_Energy_Sum=0.0;
906 double Layer2_Energy_Sum=0.0;
907 double Layer3_Energy_Sum=0.0;
908 double Layer4_Energy_Sum=0.0;
910 double Layer1_Energy_Sum2=0.0;
911 double Layer2_Energy_Sum2=0.0;
912 double Layer3_Energy_Sum2=0.0;
913 double Layer4_Energy_Sum2=0.0;
915 double Layer1_Energy_Sumpos = 0.0;
916 double Layer2_Energy_Sumpos = 0.0;
917 double Layer3_Energy_Sumpos = 0.0;
918 double Layer4_Energy_Sumpos = 0.0;
919 for(
unsigned int loc_k = 0; loc_k < associated_points.size(); loc_k++)
921 const DBCALPoint* a_point = associated_points[loc_k];
922 vector<const DBCALUnifiedHit*> associated_unifiedhits;
923 a_point->Get(associated_unifiedhits);
924 int point_size = associated_points.size();
925 int module = associated_points[loc_k]->module();
926 int layer = associated_points[loc_k]->layer();
927 int sector = associated_points[loc_k]->sector();
928 double theta = associated_points[loc_k]->theta();
929 double sintheta =
sin(theta);
930 double point_energy = associated_points[loc_k]->E();
931 int channel_per_module;
933 double z = associated_points[loc_k]->z();
934 double r = associated_points[loc_k]->r();
935 double theta_wrt_vertex = atan(r/(z-kinfitVertexZ));
936 double point_energy_sum;
937 double point_energy_over_angle = TMath::Abs(point_energy/
sin(theta_wrt_vertex));
938 point_energy_sum += point_energy_over_angle;
945 if(associated_points.size()<4) logical = 0;
950 if (layer == 1 && sector == 1) channel_per_module = 0;
951 if (layer == 1 && sector == 2) channel_per_module = 1;
952 if (layer == 1 && sector == 3) channel_per_module = 2;
953 if (layer == 1 && sector == 4) channel_per_module = 3;
954 if (layer == 2 && sector == 1) channel_per_module = 4;
955 if (layer == 2 && sector == 2) channel_per_module = 5;
956 if (layer == 2 && sector == 3) channel_per_module = 6;
957 if (layer == 2 && sector == 4) channel_per_module = 7;
958 if (layer == 3 && sector == 1) channel_per_module = 8;
959 if (layer == 3 && sector == 2) channel_per_module = 9;
960 if (layer == 3 && sector == 3) channel_per_module = 10;
961 if (layer == 3 && sector == 4) channel_per_module = 11;
962 if (layer == 4 && sector == 1) channel_per_module = 12;
963 if (layer == 4 && sector == 2) channel_per_module = 13;
964 if (layer == 4 && sector == 3) channel_per_module = 14;
965 if (layer == 4 && sector == 4) channel_per_module = 15;
966 int channel = channel_per_module + (module-1)*16;
969 if(associated_points.size()>3)
972 if(layer==1)Layer1_Energy_Sum += point_energy;
973 if(layer==2)Layer2_Energy_Sum += point_energy;
974 if(layer==3)Layer3_Energy_Sum += point_energy;
975 if(layer==4)Layer4_Energy_Sum += point_energy;
978 if(layer==1)Layer1_Energy_Sum2 += point_energy/sintheta;
979 if(layer==2)Layer2_Energy_Sum2 += point_energy/sintheta;
980 if(layer==3)Layer3_Energy_Sum2 += point_energy/sintheta;
981 if(layer==4)Layer4_Energy_Sum2 += point_energy/sintheta;
983 if(layer==1)Layer1_Energy_Sumpos += point_energy/
L1_pathlength;
984 if(layer==2)Layer2_Energy_Sumpos += point_energy/
L2_pathlength;
985 if(layer==3)Layer3_Energy_Sumpos += point_energy/
L3_pathlength;
986 if(layer==4)Layer4_Energy_Sumpos += point_energy/
L4_pathlength;
993 if(loc_k+1 == associated_points.size() ){
994 if(Layer1_Energy_Sum < Layer4_Energy_Sum){
995 if(Layer1_Energy_Sum > 0.001 && Layer2_Energy_Sum > 0.001 && Layer3_Energy_Sum > 0.001 && Layer4_Energy_Sum > 0.001)
Layer1_Energy_v2->Fill(Layer1_Energy_Sum);
996 if(Layer1_Energy_Sum > 0.001 && Layer2_Energy_Sum > 0.001 && Layer3_Energy_Sum > 0.001 && Layer4_Energy_Sum > 0.001)
Layer2_Energy_v2->Fill(Layer2_Energy_Sum);
997 if(Layer1_Energy_Sum > 0.001 && Layer2_Energy_Sum > 0.001 && Layer3_Energy_Sum > 0.001 && Layer4_Energy_Sum > 0.001)
Layer3_Energy_v2->Fill(Layer3_Energy_Sum);
998 if(Layer1_Energy_Sum > 0.001 && Layer2_Energy_Sum > 0.001 && Layer3_Energy_Sum > 0.001 && Layer4_Energy_Sum > 0.001)
Layer4_Energy_v2->Fill(Layer4_Energy_Sum);
1000 if(Layer1_Energy_Sum2 < Layer4_Energy_Sum2 && Layer1_Energy_Sum2 > 0.005){
1006 if(Layer1_Energy_Sumpos < Layer4_Energy_Sumpos && Layer1_Energy_Sumpos > 0.005 && Layer4_Energy_Sumpos > 0.005){
1017 for(
unsigned int loc_m = 0; loc_m < associated_unifiedhits.size(); loc_m++)
1019 int modulehit = associated_unifiedhits[loc_m]->module;
1020 int layerhit = associated_unifiedhits[loc_m]->layer;
1021 int sectorhit = associated_unifiedhits[loc_m]->sector;
1022 int end = associated_unifiedhits[loc_m]->end;
1023 double unifiedhit_energy = associated_unifiedhits[loc_m]->E;
1323 for(
unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
1325 const DBCALShower* locBCALShower = locBCALShowers[loc_i];
1326 if(find(matchedShowersneg.begin(), matchedShowersneg.end(), locBCALShower) == matchedShowersneg.end())
continue;
1331 double shower_x = locBCALShower->
x;
1332 double shower_y = locBCALShower->
y;
1334 r = TMath::Sqrt(TMath::Power(shower_x,2)+TMath::Power(shower_y,2));
1335 z = locBCALShower->
z;
1338 vector<const DBCALCluster*> associated_clusters;
1339 locBCALShower->Get(associated_clusters);
1341 if(associated_clusters.size() == 0) cout <<
" assoc cluster is empty " << endl;
1343 for(
unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
1346 const DBCALCluster* a_cluster = associated_clusters[loc_j];
1347 vector<const DBCALPoint*> associated_points;
1348 a_cluster->Get(associated_points);
1353 double layer1_energysum_inter = 0.0;
1354 double layer2_energysum_inter = 0.0;
1355 double layer3_energysum_inter = 0.0;
1356 double layer4_energysum_inter = 0.0;
1357 phi = associated_clusters[loc_j]->phi();
1358 theta = associated_clusters[loc_j]->theta();
1360 for(
unsigned int loc_k = 0 ; loc_k < associated_points.size(); loc_k++)
1363 module = associated_points[loc_k]->module();
1364 layer = associated_points[loc_k]->layer();
1365 sector = associated_points[loc_k]->sector();
1370 double first_pos = TMath::Sqrt( TMath::Power(myposL1.X(),2)+TMath::Power(myposL1.Y(),2)+TMath::Power(myposL1.Z(),2));
1371 double second_pos = TMath::Sqrt( TMath::Power(myposL2.X(),2)+TMath::Power(myposL2.Y(),2)+TMath::Power(myposL2.Z(),2));
1372 double third_pos = TMath::Sqrt( TMath::Power(myposL3.X(),2)+TMath::Power(myposL3.Y(),2)+TMath::Power(myposL3.Z(),2));
1373 double fourth_pos = TMath::Sqrt( TMath::Power(myposL4.X(),2)+TMath::Power(myposL4.Y(),2)+TMath::Power(myposL4.Z(),2));
1374 double fifth_pos = TMath::Sqrt( TMath::Power(myposL5.X(),2)+TMath::Power(myposL5.Y(),2)+TMath::Power(myposL5.Z(),2));
1384 if(loc_k+1 == associated_points.size()){
1392 int channel_per_module;
1393 if (
layer == 1 &&
sector == 1) channel_per_module = 0;
1394 if (
layer == 1 &&
sector == 2) channel_per_module = 1;
1395 if (
layer == 1 &&
sector == 3) channel_per_module = 2;
1396 if (
layer == 1 &&
sector == 4) channel_per_module = 3;
1397 if (
layer == 2 &&
sector == 1) channel_per_module = 4;
1398 if (
layer == 2 &&
sector == 2) channel_per_module = 5;
1399 if (
layer == 2 &&
sector == 3) channel_per_module = 6;
1400 if (
layer == 2 &&
sector == 4) channel_per_module = 7;
1401 if (
layer == 3 &&
sector == 1) channel_per_module = 8;
1402 if (
layer == 3 &&
sector == 2) channel_per_module = 9;
1403 if (
layer == 3 &&
sector == 3) channel_per_module = 10;
1404 if (
layer == 3 &&
sector == 4) channel_per_module = 11;
1405 if (
layer == 4 &&
sector == 1) channel_per_module = 12;
1406 if (
layer == 4 &&
sector == 2) channel_per_module = 13;
1407 if (
layer == 4 &&
sector == 3) channel_per_module = 14;
1408 if (
layer == 4 &&
sector == 4) channel_per_module = 15;
1422 for(
unsigned int loc_i = 0; loc_i < locBCALShowers.size(); loc_i++)
1424 const DBCALShower* locBCALShower = locBCALShowers[loc_i];
1425 if(find(matchedShowerspos.begin(), matchedShowerspos.end(), locBCALShower) == matchedShowerspos.end())
continue;
1428 double shower_x = locBCALShower->
x;
1429 double shower_y = locBCALShower->
y;
1431 z = locBCALShower->
z;
1434 r = TMath::Sqrt(TMath::Power(shower_x,2)+TMath::Power(shower_y,2));
1435 vector<const DBCALCluster*> associated_clusters;
1436 locBCALShower->Get(associated_clusters);
1438 if(associated_clusters.size() == 0) cout <<
" assoc cluster is empty " << endl;
1440 for(
unsigned int loc_j = 0; loc_j < associated_clusters.size(); loc_j++)
1443 const DBCALCluster* a_cluster = associated_clusters[loc_j];
1444 vector<const DBCALPoint*> associated_points;
1445 a_cluster->Get(associated_points);
1447 phi = associated_clusters[loc_j]->phi();
1448 theta = associated_clusters[loc_j]->theta();
1454 double layer1_energysum_inter = 0.0;
1455 double layer2_energysum_inter = 0.0;
1456 double layer3_energysum_inter = 0.0;
1457 double layer4_energysum_inter = 0.0;
1459 for(
unsigned int loc_k = 0 ; loc_k < associated_points.size(); loc_k++)
1462 module = associated_points[loc_k]->module();
1463 layer = associated_points[loc_k]->layer();
1464 sector = associated_points[loc_k]->sector();
1469 double first_pos = TMath::Sqrt( TMath::Power(myposL1.X(),2)+TMath::Power(myposL1.Y(),2)+TMath::Power(myposL1.Z(),2));
1470 double second_pos = TMath::Sqrt( TMath::Power(myposL2.X(),2)+TMath::Power(myposL2.Y(),2)+TMath::Power(myposL2.Z(),2));
1471 double third_pos = TMath::Sqrt( TMath::Power(myposL3.X(),2)+TMath::Power(myposL3.Y(),2)+TMath::Power(myposL3.Z(),2));
1472 double fourth_pos = TMath::Sqrt( TMath::Power(myposL4.X(),2)+TMath::Power(myposL4.Y(),2)+TMath::Power(myposL4.Z(),2));
1473 double fifth_pos = TMath::Sqrt( TMath::Power(myposL5.X(),2)+TMath::Power(myposL5.Y(),2)+TMath::Power(myposL5.Z(),2));
1479 int channel_per_module;
1480 if (
layer == 1 &&
sector == 1) channel_per_module = 0;
1481 if (
layer == 1 &&
sector == 2) channel_per_module = 1;
1482 if (
layer == 1 &&
sector == 3) channel_per_module = 2;
1483 if (
layer == 1 &&
sector == 4) channel_per_module = 3;
1484 if (
layer == 2 &&
sector == 1) channel_per_module = 4;
1485 if (
layer == 2 &&
sector == 2) channel_per_module = 5;
1486 if (
layer == 2 &&
sector == 3) channel_per_module = 6;
1487 if (
layer == 2 &&
sector == 4) channel_per_module = 7;
1488 if (
layer == 3 &&
sector == 1) channel_per_module = 8;
1489 if (
layer == 3 &&
sector == 2) channel_per_module = 9;
1490 if (
layer == 3 &&
sector == 3) channel_per_module = 10;
1491 if (
layer == 3 &&
sector == 4) channel_per_module = 11;
1492 if (
layer == 4 &&
sector == 1) channel_per_module = 12;
1493 if (
layer == 4 &&
sector == 2) channel_per_module = 13;
1494 if (
layer == 4 &&
sector == 3) channel_per_module = 14;
1495 if (
layer == 4 &&
sector == 4) channel_per_module = 15;
1503 if(loc_k+1 == associated_points.size()){
1525 for (
unsigned int i=0; i < locTrackTimeBased.size() ; i++)
1527 if (find(matchedTracks.begin(), matchedTracks.end(), locTrackTimeBased[i]) == matchedTracks.end())
continue;
1528 double charge = locTrackTimeBased[i]->rt->q;
1529 double p = locTrackTimeBased[i]->momentum().Mag();
1532 for(
unsigned int i=0; i<locBCALShowers.size(); i++)
1535 if (find(matchedShowers.begin(), matchedShowers.end(),s1) == matchedShowers.end())
continue;
1536 double E_MatchedShower = s1->
E;
1576 for(
size_t loc_i = 0; loc_i < locBCALShowers.size(); ++loc_i)
1578 double x = locBCALShowers[loc_i]->x;
1579 double y = locBCALShowers[loc_i]->y;
1580 double z = locBCALShowers[loc_i]->z;
1581 double E = locBCALShowers[loc_i]->E;
1583 double R = pos_bcal.Perp();
1584 double R1 =
sqrt((x-kinfitVertexX)*(x-kinfitVertexX)+(y-kinfitVertexY)*(y-kinfitVertexY)+(z-kinfitVertexZ)*(z-kinfitVertexZ));
1585 double phi = pos_bcal.Phi();
1586 for(
size_t j = 0 ; j < locTrackTimeBased.size(); ++j)
1589 locTrackTimeBased[j]->rt->GetIntersectionWithRadius(R, mypos);
1591 double dPhi = TMath::Abs(mypos.Phi()-pos_bcal.Phi());
1592 double dZ = TMath::Abs(mypos.Z() -
z);
1593 TLorentzVector
p1(E*(x-kinfitVertexX)/R, E*(y-kinfitVertexY)/R1, E*(z-kinfitVertexZ)/R1, E);
1609 for(
unsigned int j = 0 ; j < locTrackTimeBased.size(); ++j)
1611 for(
unsigned int i = 0 ; i < locFCALShowers.size(); ++i)
1613 double E = locFCALShowers[i]->getEnergy();
1614 double x = locFCALShowers[i]->getPosition().X();
1615 double y = locFCALShowers[i]->getPosition().Y();
1616 double z = locFCALShowers[i]->getPosition().Z();
1622 double myorigin = pos_fcal.Perp();
1625 locTrackTimeBased[j]->rt->GetIntersectionWithPlane(pos_fcal,norm,mypos2,mymom);
1626 double pmag = mymom.Mag();
1630 double dX = TMath::Abs(mypos2.X() -
x);
1631 double dY = TMath::Abs(mypos2.Y() -
y);
1632 double dZ = TMath::Abs(mypos2.Z() -
z);
1633 double dR = TMath::Abs(mypos2.Perp() - pos_fcal.Perp());
1648 for(
size_t j = 0 ; j < locTrackTimeBased.size(); ++j)
1650 for(
size_t loc_i = 0; loc_i < locFCALClusters.size(); ++loc_i)
1655 DVector3 fcalpos=locFCALClusters[loc_i]->getCentroid();
1663 if(locTrackTimeBased[j]->rt->GetIntersectionWithPlane(fcalpos,norm,pos,mom)==NOERROR)
1665 double diffX = TMath::Abs(fcalpos.X() - pos.X());
1666 double diffY = TMath::Abs(fcalpos.Y() - pos.Y());
1667 double track_mom = TMath::Sqrt( TMath::Power(mom.X(),2) +
1668 TMath::Power(mom.Y(),2) + TMath::Power(mom.Z(),2));
1669 double cluster_energy = locFCALClusters[loc_i]->getEnergy();
1670 double EvsP = cluster_energy/track_mom;
1671 double cluster_radius = TMath::Sqrt( TMath::Power(fcalpos.X(),2) +
1672 TMath::Power(fcalpos.Y(),2));
1676 if(diffX < 10 && diffY < 10){
1678 FCAL_Evsp->Fill(track_mom,cluster_energy);
1684 if (diffX < 5 && diffY < 5 && track_mom >1.3 && cluster_radius > 20.0 && cluster_radius < 60.0)
1804 for(
unsigned int i=0; i<locBCALShowers.size(); i++)
1808 if (find(matchedShowers.begin(), matchedShowers.end(),locBCALShowers[i]) != matchedShowers.end())
continue;
1810 vector<const DBCALCluster*> associated_clusters1;
1812 s1->Get(associated_clusters1);
1813 double shower_energy = s1->
E;
1814 double raw_shower_energy = s1->
E_raw;
1817 for(
unsigned int loc_j = 0; loc_j < associated_clusters1.size(); loc_j++)
1819 double cluster_energy = associated_clusters1[loc_j]->E();
1821 const DBCALCluster* a_cluster1 = associated_clusters1[loc_j];
1822 vector<const DBCALPoint*> associated_points1;
1823 a_cluster1->Get(associated_points1);
1824 for(
unsigned int loc_k = 0; loc_k < associated_points1.size(); loc_k++)
1826 module = associated_points1[loc_k]->module();
1827 layer = associated_points1[loc_k]->layer();
1828 sector = associated_points1[loc_k]->sector();
1829 int channel_per_module;
1830 if (
layer == 1 &&
sector == 1) channel_per_module = 0;
1831 if (
layer == 1 &&
sector == 2) channel_per_module = 1;
1832 if (
layer == 1 &&
sector == 3) channel_per_module = 2;
1833 if (
layer == 1 &&
sector == 4) channel_per_module = 3;
1834 if (
layer == 2 &&
sector == 1) channel_per_module = 4;
1835 if (
layer == 2 &&
sector == 2) channel_per_module = 5;
1836 if (
layer == 2 &&
sector == 3) channel_per_module = 6;
1837 if (
layer == 2 &&
sector == 4) channel_per_module = 7;
1838 if (
layer == 3 &&
sector == 1) channel_per_module = 8;
1839 if (
layer == 3 &&
sector == 2) channel_per_module = 9;
1840 if (
layer == 3 &&
sector == 3) channel_per_module = 10;
1841 if (
layer == 3 &&
sector == 4) channel_per_module = 11;
1842 if (
layer == 4 &&
sector == 1) channel_per_module = 12;
1843 if (
layer == 4 &&
sector == 2) channel_per_module = 13;
1844 if (
layer == 4 &&
sector == 3) channel_per_module = 14;
1845 if (
layer == 4 &&
sector == 4) channel_per_module = 15;
1848 theta = associated_points1[loc_k]->theta();
1923 for(
unsigned int i = 0; i < matchedShowers.size(); i++)
1925 double E = matchedShowers[i]->E;
1928 vector<const DBCALCluster*> associated_clusters;
1929 ms->Get(associated_clusters);
1930 for(
unsigned int j = 0 ; j < associated_clusters.size(); j++)
1932 double cluster_E = associated_clusters[j]->E();
1940 for(
unsigned int i=0; i<locBCALShowers.size(); i++)
1944 double E_before = locBCALShowers[i]->E;
1947 if (find(matchedShowers.begin(), matchedShowers.end(),locBCALShowers[i]) != matchedShowers.end())
continue;
1949 vector<const DBCALCluster*> associated_clusters1;
1950 s1->Get(associated_clusters1);
1951 double dx1 = s1->
x - kinfitVertexX;
1952 double dy1 = s1->
y - kinfitVertexY;
1953 double dz1 = s1->
z - kinfitVertexZ;
1960 double R1 =
sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
1966 double ECalc = s1->
E*(1.106+(dz1+65.0-208.4)*(dz1+65.0-208.4)*6.851E-6);
1967 TLorentzVector
p1(
E1*dx1/R1,
E1*dy1/R1,
E1*dz1/R1,
E1);
1969 double thetadeg1, thetarad1, phideg1, phirad1;
1970 thetadeg1 = p1.Theta()*180.0/TMath::Pi();
1971 thetarad1 = p1.Theta();
1972 phideg1 = p1.Phi()*180.0/TMath::Pi();
1980 TLorentzVector p1Calc(ECalc*dx1/R1, ECalc*dy1/R1, ECalc*dz1/R1, ECalc);
1983 for(
unsigned int j=i+1; j<locBCALShowers.size(); j++){
1988 if (find(matchedShowers.begin(), matchedShowers.end(),s2) != matchedShowers.end())
continue;
1990 vector<const DBCALCluster*> associated_clusters2;
1991 s2->Get(associated_clusters2);
1992 double dx2 = s2->
x - kinfitVertexX;
1993 double dy2 = s2->
y - kinfitVertexY;
1994 double dz2 = s2->
z - kinfitVertexZ;
2000 double deltaT = TMath::Abs(dt1-dt2);
2002 double R2 =
sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
2006 double ECalc = s2->
E*(1.106+(dz2+65.0-208.4)*(dz2+65.0-208.4)*6.851E-6);
2007 TLorentzVector
p2(
E2*dx2/R2,
E2*dy2/R2,
E2*dz2/R2,
E2);
2011 TLorentzVector p2Calc(ECalc*dx2/R2, ECalc*dy2/R2, ECalc*dz2/R2, ECalc);
2012 double thetadeg2, thetarad2, phideg2, phirad2, cospsi, psi;
2013 thetarad2 = p2.Theta();
2016 thetadeg2 = p2.Theta()*180.0/TMath::Pi();
2017 phideg2 = p2.Phi()*180.0/TMath::Pi();
2018 cospsi =
sin(thetarad1)*
sin(thetarad2)*cos(phirad1-phirad2)+cos(thetarad1)*cos(thetarad2);
2019 psi=acos(cospsi)*180/TMath::Pi();
2020 TLorentzVector ptot = p1+
p2;
2021 TLorentzVector ptot_raw = p1_raw + p2_raw ;
2031 double makes_sense = 0;
2032 if (R1 > R2 && dt1 > dt2) makes_sense = 1;
2033 if (R1 < R2 && dt1 < dt2) makes_sense = 1;
2034 if (R1 < R2 && dt1 > dt2) makes_sense = 0;
2035 if (R2 < R1 && dt1 > dt1) makes_sense = 0;
2090 if (
E1 > 0.5 &&
E2 > 0.5 && kinfitVertexZ > 63.0 && kinfitVertexZ < 67.0 && kinfitVertexX > -15.0 && kinfitVertexX < 15.0 && kinfitVertexY > -15.0 && kinfitVertexY < 15.0)
two_gamma_mass->Fill(ptot.M());
2091 if (
E1 > 0.2 &&
E2 > 0.2 && kinfitVertexZ > 63.0 && kinfitVertexZ < 67.0)
bcal_diphoton_mass->Fill(ptot.M());
2092 if (
E1 > 0.7 &&
E2 > 0.7 && kinfitVertexZ > 63.0 && kinfitVertexZ < 67.0 && kinfitVertexX > -15.0 && kinfitVertexX < 15.0 && kinfitVertexY > -15.0 && kinfitVertexY < 15.0)
bcal_diphoton_mass_highE->Fill(ptot.M());
2142 for(
unsigned int j=0; j<locFCALShowers.size(); j++)
2147 if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[j]) != matchedFCALShowers.end())
continue;
2150 double dx3 = s3->
getPosition().X() - kinfitVertexX;
2151 double dy3 = s3->
getPosition().Y() - kinfitVertexY;
2152 double dz3 = s3->
getPosition().Z() - kinfitVertexZ;
2154 double R3 =
sqrt(dx3*dx3 + dy3*dy3 + dz3*dz3);
2163 TLorentzVector p3(E3*dx3/R3, E3*dy3/R3, E3*dz3/R3, E3);
2165 double thetadeg3, thetarad3, phideg3, phirad3, cospsi3, psi3;
2166 thetarad3 = p3.Theta();
2168 thetadeg3 = p3.Theta()*180.0/TMath::Pi();
2169 phideg3 = p3.Phi()*180.0/TMath::Pi();
2173 for(
unsigned int i=0; i<locBCALShowers.size(); i++)
2178 if (find(matchedShowers.begin(), matchedShowers.end(),locBCALShowers[i]) != matchedShowers.end())
continue;
2180 double dx1 = s1->
x - kinfitVertexX;
2181 double dy1 = s1->
y - kinfitVertexY;
2182 double dz1 = s1->
z - kinfitVertexZ;
2184 double R1 =
sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
2196 double ECalc = s1->
E*(1.106+(dz1+65.0-208.4)*(dz1+65.0-208.4)*6.851E-6);
2197 TLorentzVector
p1(E1*dx1/R1, E1*dy1/R1, E1*dz1/R1, E1);
2199 double thetadeg1, thetarad1, phideg1, phirad1;
2200 thetadeg1 = p1.Theta()*180.0/TMath::Pi();
2201 thetarad1 = p1.Theta();
2202 phideg1 = p1.Phi()*180.0/TMath::Pi();
2205 TLorentzVector p1Calc(ECalc*dx1/R1, ECalc*dy1/R1, ECalc*dz1/R1, ECalc);
2206 TLorentzVector ptot = p1+p3;
2209 cospsi3 =
sin(thetarad1)*
sin(thetarad3)*cos(phirad1-phirad3)+cos(thetarad1)*cos(thetarad3);
2210 psi3=acos(cospsi3)*180/TMath::Pi();
2212 if( kinfitVertexZ > 63.0 && kinfitVertexZ < 67.0 && E1 > 0.5 && E3 > 0.5)
bcal_fcal_two_gamma_mass->Fill(ptot.M());
2229 for(
unsigned int j=0; j<locFCALClusters.size(); j++)
2234 if (find(matchedFCALClusters.begin(), matchedFCALClusters.end(),locFCALClusters[j]) != matchedFCALClusters.end())
continue;
2237 double dx3 = s3->
getCentroid().X() - kinfitVertexX;
2238 double dy3 = s3->
getCentroid().Y() - kinfitVertexY;
2239 double dz3 = s3->
getCentroid().Z() - kinfitVertexZ;
2241 double R3 =
sqrt(dx3*dx3 + dy3*dy3 + dz3*dz3);
2250 TLorentzVector p3(E3*dx3/R3, E3*dy3/R3, E3*dz3/R3, E3);
2252 double thetadeg3, thetarad3, phideg3, phirad3, cospsi3, psi3;
2253 thetarad3 = p3.Theta();
2255 thetadeg3 = p3.Theta()*180.0/TMath::Pi();
2256 phideg3 = p3.Phi()*180.0/TMath::Pi();
2260 for(
unsigned int i=0; i<locBCALShowers.size(); i++)
2265 if (find(matchedShowers.begin(), matchedShowers.end(),locBCALShowers[i]) != matchedShowers.end())
continue;
2267 double dx1 = s1->
x - kinfitVertexX;
2268 double dy1 = s1->
y - kinfitVertexY;
2269 double dz1 = s1->
z - kinfitVertexZ;
2271 double R1 =
sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
2283 double ECalc = s1->
E*(1.106+(dz1+65.0-208.4)*(dz1+65.0-208.4)*6.851E-6);
2284 TLorentzVector
p1(E1_raw*dx1/R1, E1_raw*dy1/R1, E1_raw*dz1/R1, E1_raw);
2286 double thetadeg1, thetarad1, phideg1, phirad1;
2287 thetadeg1 = p1.Theta()*180.0/TMath::Pi();
2288 thetarad1 = p1.Theta();
2289 phideg1 = p1.Phi()*180.0/TMath::Pi();
2292 TLorentzVector p1Calc(ECalc*dx1/R1, ECalc*dy1/R1, ECalc*dz1/R1, ECalc);
2293 TLorentzVector ptot = p1+p3;
2296 cospsi3 =
sin(thetarad1)*
sin(thetarad3)*cos(phirad1-phirad3)+cos(thetarad1)*cos(thetarad3);
2297 psi3=acos(cospsi3)*180/TMath::Pi();
2319 for(
unsigned int ii=0; ii<locFCALShowers.size(); ii++)
2325 if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[ii]) != matchedFCALShowers.end())
continue;
2329 double dx1 = s2->
getPosition().X() - kinfitVertexX;
2330 double dy1 = s2->
getPosition().Y() - kinfitVertexY;
2331 double dz1 = s2->
getPosition().Z() - kinfitVertexZ;
2332 double R1 =
sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
2336 TLorentzVector
p1(
E1*dx1/R1,
E1*dy1/R1,
E1*dz1/R1,
E1);
2348 double thetarad2 = p1.Theta();
2349 double phirad2 = p1.Phi();
2350 double thetadeg2 = p1.Theta()*180.0/TMath::Pi();
2351 double phideg2 = p1.Phi()*180.0/TMath::Pi();
2357 for(
unsigned int k=ii+1; k<locFCALShowers.size(); k++)
2361 if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[k]) != matchedFCALShowers.end())
continue;
2363 double dx2 = s3->
getPosition().X() - kinfitVertexX;
2364 double dy2 = s3->
getPosition().Y() - kinfitVertexY;
2366 double R2 =
sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
2375 TLorentzVector
p2(
E2*dx2/R2,
E2*dy2/R2,
E2*dz2/R2,
E2);
2377 double thetadeg3,thetarad3,phideg3,phirad3,cospsi3,psi3;
2378 int makes_sense3 = 0;
2379 if(R2 > R1 && dt2 > dt1) makes_sense3=1;
2380 if(R2 < R1 && dt2 < dt1) makes_sense3=1;
2381 if(R2 > R1 && dt2 < dt1) makes_sense3=0;
2382 if(R2 < R1 && dt2 > dt1) makes_sense3=0;
2383 TLorentzVector ptot = p2+
p1;
2385 thetadeg3 = p2.Theta()*180.0/TMath::Pi();
2386 thetarad3 = p2.Theta();
2387 phideg3 = p2.Phi()*180.0/TMath::Pi();
2389 cospsi3 =
sin(thetarad2)*
sin(thetarad3)*cos(phirad2-phirad3)+cos(thetarad2)*cos(thetarad3);
2390 psi3=acos(cospsi3)*180/TMath::Pi();
2401 if( kinfitVertexZ > 62.0 && kinfitVertexZ < 68.0 && E2 > 0.5 &&
E1 > 0.5 && kinfitVertexX < 15.0 && kinfitVertexX > -15.0 && kinfitVertexY > -15.0 && kinfitVertexY < 15.0)
two_fcal_gamma_mass->Fill(ptot.M());
2409 for(
unsigned int ii=0; ii<locFCALShowers.size(); ii++)
2415 if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[ii]) != matchedFCALShowers.end())
continue;
2419 double dx1 = s2->
getPosition().X() - kinfitVertexX;
2420 double dy1 = s2->
getPosition().Y() - kinfitVertexY;
2421 double dz1 = s2->
getPosition().Z() - kinfitVertexZ;
2422 double R1 =
sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
2426 TLorentzVector
p1(
E1*dx1/R1,
E1*dy1/R1,
E1*dz1/R1,
E1);
2438 double thetarad2 = p1.Theta();
2439 double phirad2 = p1.Phi();
2440 double thetadeg2 = p1.Theta()*180.0/TMath::Pi();
2441 double phideg2 = p1.Phi()*180.0/TMath::Pi();
2447 for(
unsigned int k=ii+1; k<locFCALShowers.size(); k++)
2451 if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[k]) != matchedFCALShowers.end())
continue;
2453 double dx2 = s3->
getPosition().X() - kinfitVertexX;
2454 double dy2 = s3->
getPosition().Y() - kinfitVertexY;
2455 double dz2 = s3->
getPosition().Z() - kinfitVertexZ;
2456 double R2 =
sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
2465 TLorentzVector
p2(
E2*dx2/R2,
E2*dy2/R2,
E2*dz2/R2,
E2);
2467 double thetadeg3,thetarad3,phideg3,phirad3,cospsi3,psi3;
2468 int makes_sense3 = 0;
2469 if(R2 > R1 && dt2 > dt1) makes_sense3=1;
2470 if(R2 < R1 && dt2 < dt1) makes_sense3=1;
2471 if(R2 > R1 && dt2 < dt1) makes_sense3=0;
2472 if(R2 < R1 && dt2 > dt1) makes_sense3=0;
2473 thetadeg3 = p2.Theta()*180.0/TMath::Pi();
2474 thetarad3 = p2.Theta();
2475 phideg3 = p2.Phi()*180.0/TMath::Pi();
2477 cospsi3 =
sin(thetarad2)*
sin(thetarad3)*cos(phirad2-phirad3)+cos(thetarad2)*cos(thetarad3);
2478 psi3=acos(cospsi3)*180/TMath::Pi();
2484 for(
unsigned int m = ii+2; m < locFCALShowers.size(); m ++)
2486 if (find(matchedFCALShowers.begin(), matchedFCALShowers.end(),locFCALShowers[m]) != matchedFCALShowers.end())
continue;
2487 double dx3 = locFCALShowers[m]->getPosition().X() - kinfitVertexX;
2488 double dy3 = locFCALShowers[m]->getPosition().Y() - kinfitVertexY;
2489 double dz3 = locFCALShowers[m]->getPosition().Z() - kinfitVertexZ;
2490 double R3 =
sqrt(dx3*dx3+dy3*dy3+dz3*dz3);
2491 E3 = locFCALShowers[m]->getEnergy();
2492 double dt3 = locFCALShowers[m]->getTime();
2493 x3 = locFCALShowers[m]->getPosition().X();
2494 y3 = locFCALShowers[m]->getPosition().Y();
2495 z3 = locFCALShowers[m]->getPosition().Z();
2496 t3 = locFCALShowers[m]->getTime();
2497 TLorentzVector p3(
E3*dx3/R3,
E3*dy3/R3,
E3*dz3/R3,
E3);
2499 TLorentzVector ptot = p1 + p2 + p3;
2514 for(
unsigned int ii=0; ii<locFCALClusters.size(); ii++)
2520 if (find(matchedFCALClusters.begin(), matchedFCALClusters.end(),locFCALClusters[ii]) != matchedFCALClusters.end())
continue;
2524 double dx1 = s2->
getCentroid().X() - kinfitVertexX;
2525 double dy1 = s2->
getCentroid().Y() - kinfitVertexY;
2526 double dz1 = s2->
getCentroid().Z() - kinfitVertexZ;
2527 double R1 =
sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
2531 TLorentzVector
p1(
E1*dx1/R1,
E1*dy1/R1,
E1*dz1/R1,
E1);
2543 double thetarad2 = p1.Theta();
2544 double phirad2 = p1.Phi();
2545 double thetadeg2 = p1.Theta()*180.0/TMath::Pi();
2546 double phideg2 = p1.Phi()*180.0/TMath::Pi();
2552 for(
unsigned int k=ii+1; k<locFCALClusters.size(); k++)
2556 if (find(matchedFCALClusters.begin(), matchedFCALClusters.end(),locFCALClusters[k]) != matchedFCALClusters.end())
continue;
2558 double dx2 = s3->
getCentroid().X() - kinfitVertexX;
2559 double dy2 = s3->
getCentroid().Y() - kinfitVertexY;
2560 double dz2 = s3->
getCentroid().Z() - kinfitVertexZ;
2561 double R2 =
sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
2570 TLorentzVector
p2(
E2*dx2/R2,
E2*dy2/R2,
E2*dz2/R2,
E2);
2572 double thetadeg3,thetarad3,phideg3,phirad3,cospsi3,psi3;
2573 int makes_sense3 = 0;
2574 if(R2 > R1 && dt2 > dt1) makes_sense3=1;
2575 if(R2 < R1 && dt2 < dt1) makes_sense3=1;
2576 if(R2 > R1 && dt2 < dt1) makes_sense3=0;
2577 if(R2 < R1 && dt2 > dt1) makes_sense3=0;
2578 thetadeg3 = p2.Theta()*180.0/TMath::Pi();
2579 thetarad3 = p2.Theta();
2580 phideg3 = p2.Phi()*180.0/TMath::Pi();
2582 cospsi3 =
sin(thetarad2)*
sin(thetarad3)*cos(phirad2-phirad3)+cos(thetarad2)*cos(thetarad3);
2583 psi3=acos(cospsi3)*180/TMath::Pi();
2589 TLorentzVector ptot = p1 +
p2;
2702 double Ehit_tot = 0.0;
2703 for(
unsigned int i=0; i<bcalhits.size(); i++){
2704 Ehit_tot += bcalhits[i]->E;
2709 double Etruth_tot = 0.0;
2710 double z_truth = 0.0;
2711 for(
unsigned int i=0; i<truthshowers.size(); i++){
2712 Etruth_tot += truthshowers[i]->E;
2713 z_truth += truthshowers[i]->E*truthshowers[i]->z;
2715 z_truth/=Etruth_tot;
2719 double Etot_thrown=0.0;
2720 for(
unsigned int i=0; i<mcthrowns.size(); i++){
2721 Etot_thrown += mcthrowns[i]->energy();
2722 for(
unsigned int j=0; j<locBCALShowers.size(); j++){
2723 double z = locBCALShowers[j]->z;
2726 double E = locBCALShowers[j]->E*(1.106+(z-208.4)*(z-208.4)*6.851E-6);
2735 if(mcthrowns.size()==1){
2736 const DMCThrown* mcthrown = mcthrowns[0];
2737 if(mcthrown->
momentum().Theta()>0.0001){
2739 double Ethrown = 1.0;
TH2F * Layer1_Energy_vs_Channel
DVector3 getCentroid() const
const DAnalysisUtilities * dAnalysisUtilities
jerror_t init(void)
Called once at program start.
uint32_t FCALClusters_per_event
TTree * Split_Gamma_Neutrals_raw
TH1F * Eoverp_minus_nocuts
const DEventWriterREST * dEventWriterREST
TH1F * Phi_Hist_Both_BCAL
jerror_t brun(jana::JEventLoop *locEventLoop, int locRunNumber)
Called every time a new run number is detected.
TH1F * Theta_Hist_Both_FCAL
TTree * BCALPoint_Neutral
TTree * FCALClusterNeutrals
uint32_t FCALShowers_per_event
TH1F * Psi_Hist_Both_BCAL
const DVector3 & position(void) const
TH1F * two_fcal_gamma_mass
uint32_t BCALShowers_per_event
jerror_t erun(void)
Called every time run number changes, provided brun has been called.
TH1F * Phi_Hist_Both_FCAL
TTree * BCALPoint_Charged_pos
TH1F * bcal_diphoton_mass
double Calc_DOCAVertex(const DVector3 &locUnitDir1, const DVector3 &locUnitDir2, const DVector3 &locVertex1, const DVector3 &locVertex2, DVector3 &locDOCAPoint) const
TH1F * Psi_Hist_Split_Gammas
TH1F * Theta_Hist_Both_BCAL
TH2F * Layer4_Energy_vs_Channel
TTree * BCALPoint_Charged_neg
TH2F * Layer2_Energy_vs_Channel
TH1F * Eoverp_plus_nocuts
TH1F * FCAL_Eoverp_nocuts
jerror_t fini(void)
Called after last event of last event source has been processed.
TH1F * Phi_Hist_Split_Gammas
TH1F * Thrown_Vertex_Frequency
Float_t energy_raw_shower
const DVector3 & momentum(void) const
TTree * Split_Gamma_Neutrals
TH1F * Psi_Hist_Both_FCAL
TH2F * Layer3_Energy_vs_Channel
TH1F * bcal_fcal_two_gamma_mass
TH1F * bcal_diphoton_mass_highE
TH1F * BCALShowerTrack_Energy
DVector3 getPosition() const
jerror_t evnt(jana::JEventLoop *locEventLoop, uint64_t locEventNumber)
Called every event.
TH1F * two_fcal_gamma_mass_test
TTree * Triple_FCAL_Neutrals
TH1F * Thrown_Gamma_Theta
const DEventWriterROOT * dEventWriterROOT
TH1F * Theta_Hist_Split_Gammas