Partial commit with key changes as follows:
- wireOffset is now 4 not 3. Phew. Gosh Golly. Change in Armory/ClassPW.h - Concomitant changes made to Armory/PC_Stepladder_Correction.h. In this particular version, a shift of 2 anodes in run*.sh has been replaced by an offset in z - The above needs to be worked out carefully. - The rho_vs_z fixes have NOT been corrected for. This version will likely yield nonsensical results, hence. To be fixed by a subsequent push that also cleans up the histogram plotting somewhat. *Read again* The sort uses rho from the 3 wireOffset cases in MakeVertex.C The code assumes wireOffset=4. Anywhere this matters, this version will break - caveat Emptor. Pushing here to plant a flag, but to also make sure the key changes are communicated across branches.
This commit is contained in:
parent
1e0af0fe9d
commit
7be45f35df
|
|
@ -64,7 +64,7 @@ public:
|
|||
|
||||
Coord Crossover[24][24][2];
|
||||
|
||||
|
||||
inline TVector3 getClosestWirePosAtWirePhi(std::pair<TVector3, TVector3>, double phi);
|
||||
inline std::tuple<std::pair<TVector3, TVector3>, double, double, double> GetPseudoWire(const std::vector<std::tuple<int,double,double>>& cluster, std::string type);
|
||||
|
||||
inline std::tuple<TVector3,double,double,double,double,double,double,double>
|
||||
|
|
@ -116,7 +116,7 @@ private:
|
|||
TVector3 trackVec;
|
||||
|
||||
const int nWire = 24;
|
||||
const int wireShift = 3;
|
||||
const int wireShift = 4;
|
||||
//const float zLen = 380; // mm
|
||||
// const float zLen = 348.6; // mm
|
||||
const float zLen = 174.3*2; // mm
|
||||
|
|
@ -229,6 +229,41 @@ inline void PW::ConstructGeo()
|
|||
cathodeLength = TMath::Sqrt(zLen * zLen + TMath::Power(2 * radiusC * TMath::Sin(dAngle / 2), 2)); //chord length subtending an angle alpha is 2rsin(alpha/2)
|
||||
}
|
||||
|
||||
inline TVector3 PW::getClosestWirePosAtWirePhi(std::pair<TVector3, TVector3> awire, double sx3phi_radian) {
|
||||
// 1. Get wire geometry
|
||||
TVector3 a1 = awire.first; // Top of the wire
|
||||
TVector3 a2 = awire.second; // Bottom of the wire
|
||||
TVector3 wireVec = a2 - a1; // Vector pointing down the wire
|
||||
|
||||
// Variables to track our minimums during the scan
|
||||
double min_delta_phi = 9999.0;
|
||||
double best_t = -1.0;
|
||||
TVector3 best_pcz_intersect;
|
||||
|
||||
// 2. THE SCAN: Walk down the wire in 1000 tiny steps
|
||||
// (For a 380mm wire, this is checking every 0.38 mm)
|
||||
int num_steps = 1000;
|
||||
for (int i = 0; i <= num_steps; ++i)
|
||||
{
|
||||
double t_test = (double)i / num_steps; // Ranges from 0.0 to 1.0
|
||||
TVector3 test_pt = a1 + t_test * wireVec; // The 3D point at this step
|
||||
|
||||
// Calculate absolute Delta Phi between Si hit and this specific point on the wire
|
||||
if(TMath::IsNaN(sx3phi_radian-test_pt.Phi())) continue;
|
||||
double dPhi = TMath::Abs(TVector2::Phi_mpi_pi(sx3phi_radian - test_pt.Phi())); //Phi_mpi_pi just puts the angle in the range -180 to 180
|
||||
|
||||
// If this is the smallest Delta Phi we've seen so far, save it!
|
||||
if (dPhi < min_delta_phi)
|
||||
{
|
||||
min_delta_phi = dPhi;
|
||||
best_t = t_test;
|
||||
best_pcz_intersect = test_pt;
|
||||
}
|
||||
}
|
||||
return best_pcz_intersect;
|
||||
}
|
||||
|
||||
|
||||
inline std::vector<std::vector<std::tuple<int,double,double>>>
|
||||
PW::Make_Clusters(std::unordered_map<int,std::tuple<int,double,double>> wireEvents) {
|
||||
std::vector<std::vector<std::tuple<int,double,double>>> wireClusters;
|
||||
|
|
|
|||
|
|
@ -1,36 +1,83 @@
|
|||
#include <TF1.h>
|
||||
double model(double* x, double* p) {
|
||||
/*double model(double* x, double* p) {
|
||||
double result = x[0];
|
||||
double factor = 29.0;
|
||||
double slope = 0.7;
|
||||
if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope;
|
||||
else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor;
|
||||
else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2;
|
||||
else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2;
|
||||
else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*3;
|
||||
return result;
|
||||
}
|
||||
|
||||
double model_invert(double *y, double *q) {
|
||||
double result=y[0];
|
||||
double slope = 0.7;
|
||||
double factor = 40.0;
|
||||
double factor = 0.0;
|
||||
if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope;
|
||||
else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor;
|
||||
else if(TMath::Abs(y[0]) < 85.2/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2;
|
||||
else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2;
|
||||
else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*3;
|
||||
return result;
|
||||
}*/
|
||||
|
||||
double model_invert(double* y, double* p) {
|
||||
double result = y[0];
|
||||
double slope = 0.6;
|
||||
double z_grid[8] = {147.998,101.946,59.7634,19.6965,-19.6965,-59.7634,-101.946,-147.998};
|
||||
for(int i=0;i<7;i++) {
|
||||
if(y[0] <= z_grid[i] && y[0] > z_grid[i+1]) {
|
||||
double zavg = (z_grid[i] + z_grid[i+1])*0.5; //midpoint about which we pivot
|
||||
result = (y[0]-zavg)/slope + zavg;
|
||||
break;
|
||||
}
|
||||
}
|
||||
return result+80;
|
||||
}
|
||||
|
||||
double model_a1c1(double* x, double* p) {
|
||||
double result = x[0];
|
||||
double factor = 29.0;
|
||||
double slope = 0.0;
|
||||
if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope;
|
||||
else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor;
|
||||
else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2;
|
||||
else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*3;
|
||||
return result;
|
||||
}
|
||||
|
||||
double model_invert_a1c1(double *y, double *q) {
|
||||
double result=y[0];
|
||||
/* double slope = 1.0;
|
||||
double factor = 5.0;
|
||||
if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope;
|
||||
else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor;
|
||||
else if(TMath::Abs(y[0]) < 85.2/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2;
|
||||
else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor**/;
|
||||
return result+40;
|
||||
}
|
||||
|
||||
|
||||
/*void testmodel() {
|
||||
TF1 eqline("x","x",-200,200);
|
||||
eqline.Draw("");
|
||||
eqline.SetLineStyle(kDashed);
|
||||
|
||||
//TF1 f1("model",model,-200,200,2);
|
||||
TF1 f1("model_inv",model_invert,-200,200,2);
|
||||
//TF1 f1("model",model,-200,200,2);
|
||||
TF1 f1a("model_inv",model_a1c1,-200,200,2);
|
||||
eqline.SetNpx(10000);
|
||||
f1a.SetNpx(10000);
|
||||
std::vector<double> pars = {0.0,1.};
|
||||
f1a.SetParameters(pars.data());
|
||||
f1a.SetLineColor(kGreen+2);
|
||||
f1a.SetLineStyle(kLine);
|
||||
f1a.Draw("L SAME");
|
||||
|
||||
TF1 f1("model",model,-200,200,2);
|
||||
//TF1 f1("model_inv",model_invert,-200,200,2);
|
||||
eqline.SetNpx(10000);
|
||||
f1.SetNpx(10000);
|
||||
std::vector<double> pars = {0.0,1.};
|
||||
//std::vector<double> pars = {0.0,1.};
|
||||
f1.SetParameters(pars.data());
|
||||
f1.SetLineColor(kGreen+2);
|
||||
f1.SetLineStyle(kLine);
|
||||
|
|
@ -38,5 +85,4 @@ double model_invert(double *y, double *q) {
|
|||
|
||||
gPad->Modified(); gPad->Update();
|
||||
while(gPad->WaitPrimitive());
|
||||
|
||||
}*/
|
||||
|
|
|
|||
153
MakeVertex.C
153
MakeVertex.C
|
|
@ -45,6 +45,7 @@ const double anode_gain = 1.5146e-5; //channels --> MeV
|
|||
std::string dataset;
|
||||
|
||||
TF1 pcfix_func("func",model_invert,-200,200);
|
||||
TF1 pcfix_func_a1c1("func_a1c1",model_invert_a1c1,-200,200);
|
||||
TGraph *MeV_to_cm=NULL,*cm_to_MeV=NULL;
|
||||
TGraph *MeV_to_cm_p=NULL,*cm_to_MeVp=NULL;
|
||||
|
||||
|
|
@ -411,9 +412,12 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
double backE = det.backE*sx3BackGain[id][det.stripF][det.stripB];
|
||||
//if(backE<2000) continue;
|
||||
det.stripF=3-det.stripF;
|
||||
double beta_n = 15.0 + TMath::ATan2((2*det.stripF-3)*40.30, 8.0*88.0*TMath::Cos(15.0*M_PI/180.0))*180./M_PI; //how much to add per strip to the starting position
|
||||
|
||||
double alpha_n = TMath::ATan2((2*det.stripF-3)*40.30, 8.0*88.0*TMath::Cos(15.0*M_PI/180.0))*180./M_PI; //angle subtended w.r.t the radial perpendicular bisector of each sx3
|
||||
double beta_n = 15.0+alpha_n; //how much to add per strip to the starting position? this is the angle w.r.t an edge of the sx3, the above values run as (-10.08deg, -3.39deg, 3.39deg, 10.08deg)
|
||||
double phi_n = ((-id+0.5)*30+beta_n);
|
||||
phi_n+=45;
|
||||
double rho_at_strip = 88.0/TMath::Cos(alpha_n*M_PI/180.0); //*TMath::Cos(15.0*M_PI/180.0) if the edge-length is 88mm
|
||||
|
||||
//if(getenv("flip180")) {
|
||||
// if(std::string(getenv("flip180"))=="1") {
|
||||
|
|
@ -422,7 +426,8 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
//}
|
||||
//}
|
||||
phi_n*=M_PI/180.; //starting-position phi + strip contribution
|
||||
Event sx3ev(TVector3(88.0*TMath::Cos(phi_n),88.0*TMath::Sin(phi_n),z),backE*0.001,-1,det.ts,-1,det.stripB+4*id,det.stripF+4*id);
|
||||
//Event sx3ev(TVector3(88.0*TMath::Cos(phi_n),88.0*TMath::Sin(phi_n),z),backE*0.001,-1,det.ts,-1,det.stripB+4*id,det.stripF+4*id);
|
||||
Event sx3ev(TVector3(rho_at_strip*TMath::Cos(phi_n),rho_at_strip*TMath::Sin(phi_n),z),backE*0.001,-1,det.ts,-1,det.stripB+4*id,det.stripF+4*id);
|
||||
SX3_Events.push_back(sx3ev);
|
||||
plotter->Fill2D("sx3backs_gm",100,0,100,800,0,8192,det.stripB+4*id,backE);
|
||||
|
||||
|
|
@ -442,7 +447,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
int qqqCount = 0;
|
||||
int qqqAdjCh = 0;
|
||||
// REMOVE WHEN RERUNNING USING THE NEW CALIBRATION FILE
|
||||
std::vector<Event> QQQ_Events, PC_Events;
|
||||
std::vector<Event> QQQ_Events, PC_Events,PC_Events_OnlyAnode, PC_Events_OnlyCathode;
|
||||
std::vector<Event> QQQ_Events_Raw, PC_Events_Raw;
|
||||
std::vector<Event> QQQ_Events2; //clustering done
|
||||
|
||||
|
|
@ -739,11 +744,9 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
std::vector<std::pair<double,double>> sumE_AC;
|
||||
for(auto aCluster: aClusters) {
|
||||
for(auto cCluster: cClusters) {
|
||||
if(aCluster.size() ==0 ) continue;
|
||||
if(cCluster.size() ==0 ) continue;
|
||||
//both have at least 1, here. Keep the a1, c1 events
|
||||
auto [crossover,alpha,apSumE,cpSumE,apMaxE,cpMaxE,apTSMaxE,cpTSMaxE] = pwinstance.FindCrossoverProperties(aCluster, cCluster);
|
||||
if(alpha!=9999999 && apSumE!=-1) {
|
||||
if(alpha!=9999999 && apSumE!=-1 && aCluster.size() && cCluster.size()) {
|
||||
//Event PCEvent(crossover,apMaxE,cpMaxE,apTSMaxE,cpTSMaxE);
|
||||
//Event PCEvent(crossover,apSumE,cpSumE,apTSMaxE,cpTSMaxE);
|
||||
Event PCEvent(crossover,apSumE,cpMaxE,apTSMaxE,cpTSMaxE); //run12 shows cathode-max and anode-sum provide best dE signals.
|
||||
|
|
@ -752,7 +755,21 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
PCEvent.multi2=cCluster.size();
|
||||
PC_Events.push_back(PCEvent);
|
||||
sumE_AC.push_back(std::pair(apSumE,cpSumE));
|
||||
} else {
|
||||
} else if(alpha!=9999999){
|
||||
std::cout << aCluster.size() << " " << cCluster.size() << std::endl;
|
||||
//Store anode-only events in a separate structure for later analysis
|
||||
if(aCluster.size()!=0) {
|
||||
Event PCEvent_OnlyA(crossover,apSumE,cpMaxE,apTSMaxE,cpTSMaxE); //run12 shows cathode-max and anode-sum provide best dE signals.
|
||||
PCEvent_OnlyA.multi1=aCluster.size();
|
||||
PCEvent_OnlyA.multi2=0;
|
||||
PC_Events_OnlyAnode.push_back(PCEvent_OnlyA);
|
||||
}
|
||||
if(cCluster.size()!=0) {
|
||||
Event PCEvent_OnlyC(crossover,apSumE,cpMaxE,apTSMaxE,cpTSMaxE); //run12 shows cathode-max and anode-sum provide best dE signals.
|
||||
PCEvent_OnlyC.multi1=0;
|
||||
PCEvent_OnlyC.multi2=cCluster.size();
|
||||
PC_Events_OnlyCathode.push_back(PCEvent_OnlyC);
|
||||
}
|
||||
;//std::cout << "AAAA " << std::endl;
|
||||
}
|
||||
}
|
||||
|
|
@ -781,6 +798,35 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
}
|
||||
|
||||
for(auto sx3event: SX3_Events) {
|
||||
|
||||
for(const auto acluster: aClusters) {
|
||||
auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(acluster,"ANODE");
|
||||
int a_number = acluster.size();
|
||||
TVector3 pc_closest = pwinstance.getClosestWirePosAtWirePhi(apwire,sx3event.pos.Phi());
|
||||
plotter->Fill1D("dt_anode_interp_sx3",800,-2000,2000,sx3event.Time1 - apTSMaxE,"ainterp_noc");
|
||||
|
||||
if(sx3event.Time1 - apTSMaxE < -150) {
|
||||
plotter->Fill1D("dt_anode_ainterp_sx3_gated",800,-2000,2000,sx3event.Time1 - apTSMaxE,"ainterp_noc");
|
||||
plotter->Fill2D("dEa_ainterp_Esx3_TC1_ignC"+std::to_string(acluster.size()),400,0,10,800,0,40000,sx3event.Energy1,apSumE,"ainterp_noc");
|
||||
plotter->Fill2D("pcPhi_ainterp_sx3Phi_TC1_ignC"+std::to_string(acluster.size()),120,-360,360,120,-360,360,pc_closest.Phi()*180./M_PI,sx3event.pos.Phi()*180./M_PI,"ainterp_noc");
|
||||
plotter->Fill2D("pcZ_ainterp_sx3Z_TC1_ignC"+std::to_string(acluster.size()),300,-100,200,400,-200,200,sx3event.pos.Z(),pc_closest.Z(),"ainterp_noc");
|
||||
|
||||
double sx3theta = TMath::ATan2(88.0,sx3event.pos.Z()-source_vertex);
|
||||
double pczguess = 37.0/TMath::Tan(sx3theta) + source_vertex;
|
||||
double sinTheta = TMath::Sin(sx3theta);
|
||||
plotter->Fill2D("pcZ_ainterp_sx3pczguess_TC1_ignC"+std::to_string(acluster.size()),300,-100,200,400,-200,200,pczguess,pc_closest.Z(),"ainterp_noc");
|
||||
|
||||
|
||||
TVector3 x2(pc_closest), x1(sx3event.pos);
|
||||
TVector3 v = x2-x1;
|
||||
double t_minimum = -1.0*(x1.X()*v.X()+x1.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y());
|
||||
TVector3 vector_closest_to_axis = x1 + t_minimum*v;
|
||||
plotter->Fill2D("vertexZ_ainterp_sx3Z_TC1_ignC"+std::to_string(acluster.size()),300,-100,200,400,-200,200,sx3event.pos.Z(),vector_closest_to_axis.Z(),"ainterp_noc");
|
||||
plotter->Fill2D("vertexXY_ainterp_TC1_ignC"+std::to_string(acluster.size()),200,-100,100,200,-100,100,vector_closest_to_axis.X(),vector_closest_to_axis.Y(),"ainterp_noc");
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
for(int i=0; i<24; i++) {
|
||||
if(aWireEvents.find(i) != aWireEvents.end()) {
|
||||
auto awire = aWireEvents[i];
|
||||
|
|
@ -789,6 +835,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
//plotter->Fill2D("sx3_z_strip#_awire"+std::to_string(std::get<0>(awire)), 400,-100,100, 100, -50,50,sx3event.pos.Z(), sx3event.ch2);
|
||||
plotter->Fill2D("onewire_dEa_Esx3_TC1_fullev"+std::to_string(PC_Events.size()>0),400,0,10,800,0,40000,sx3event.Energy1,std::get<1>(awire));
|
||||
plotter->Fill2D("onewire_aNum_sx3Phi_TC1_fullev"+std::to_string(PC_Events.size()>0),24,0,24,120,-360,360,i,sx3event.pos.Phi()*180./M_PI);
|
||||
//plotter->Fill2D("sx3_z_phi_ow_awire"+std::to_string(anodeIndex)+"_sx3strip"+std::to_string(sx3event.ch2), 400,-100,100, 200, -200,200,sx3event.pos.Z(), sx3event.pos.Phi()*180/M_PI );
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -824,18 +871,9 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
}//for 'i' loop
|
||||
}
|
||||
|
||||
|
||||
for(auto pcevent:PC_Events) {
|
||||
if(aClusters.size()==1 && cClusters.size() == 1) {
|
||||
//plotter->Fill1D("pcz_a"+std::to_string(aClusters.at(0).size())+"_c"+std::to_string(cClusters.at(0).size()),800,-200,200,pcevent.pos.Z(),"wiremult");
|
||||
std::string detid="_+_";
|
||||
if(SX3_Events.size()) detid="+sx3";
|
||||
if(QQQ_Events.size()) detid="+qqq";
|
||||
//plotter->Fill1D("pcz_a"+std::to_string(aClusters.at(0).size())+"_c"+std::to_string(cClusters.at(0).size())+detid,800,-200,200,pcevent.pos.Z(),"wiremult");
|
||||
}
|
||||
|
||||
for(auto sx3event:SX3_Events) {
|
||||
PCSX3TimeCut=false;
|
||||
for(auto sx3event:SX3_Events) {
|
||||
for(auto pcevent:PC_Events) {
|
||||
plotter->Fill1D("dt_pcA_sx3B"+std::to_string(sx3event.ch2),640,-2000,2000,sx3event.Time1 - pcevent.Time1,"hTiming");
|
||||
plotter->Fill1D("dt_pcC_sx3B"+std::to_string(sx3event.ch2),640,-2000,2000,sx3event.Time1 - pcevent.Time2,"hTiming");
|
||||
if(sx3event.Time1 - pcevent.Time1 < 0)//-150 for alphas
|
||||
|
|
@ -876,6 +914,8 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
double pczguess = 37.0/TMath::Tan(sx3theta) + source_vertex;
|
||||
double pcz_guess_int = z_to_crossover_rho(pcevent.pos.Z())/TMath::Tan(sx3theta) + source_vertex;
|
||||
double sinTheta = TMath::Sin(sx3theta);
|
||||
plotter->Fill2D("dE2_E_Anodesx3B",400,0,30,800,0,40000,sx3event.Energy1,pcevent.Energy1*TMath::Sin(sx3theta));
|
||||
plotter->Fill2D("dE2_E_Cathodesx3B",400,0,30,800,0,10000,sx3event.Energy1,pcevent.Energy2*TMath::Sin(sx3theta));
|
||||
|
||||
TVector3 x2(pcevent.pos), x1(sx3event.pos);
|
||||
TVector3 v = x2-x1;
|
||||
|
|
@ -891,8 +931,25 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
|
||||
|
||||
plotter->Fill2D("pcz_vs_sx3pczguess",600,-200,200,600,-200,200,pczguess,pcevent.pos.Z()); //x-axis is all Si det, y-axis is PC anode+cathode only
|
||||
if(pcevent.multi1 == 1 && pcevent.multi2==1) {
|
||||
plotter->Fill2D("pcz_vs_sx3pczguess_A1C1",600,-200,200,600,-200,200,pczguess,pcevent.pos.Z());
|
||||
double pcz_fix_a1c1 = pcfix_func_a1c1.Eval(pcevent.pos.Z());
|
||||
|
||||
TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix_a1c1);
|
||||
TVector3 v = x2f-x1;
|
||||
double t_minimum = -1.0*(x1.X()*v.X()+x1.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y());
|
||||
TVector3 r_rhoMin_fix = x1 + t_minimum*v;
|
||||
plotter->Fill1D("VertexRecon_pczfix_sx3_a1c1",800,-300,300,r_rhoMin_fix.Z());
|
||||
plotter->Fill1D("pczfix_A1C1_1d_sx3",600,-200,200,pcz_fix_a1c1);
|
||||
plotter->Fill2D("pczfix_vs_sx3pczguess_A1C1",600,-200,200,600,-200,200,pczguess,pcz_fix_a1c1);
|
||||
plotter->Fill2D("pcz_vs_sx3pczguess_A1C1_strip"+std::to_string(sx3event.ch2),300,-200,200,600,-200,200,pczguess,pcevent.pos.Z());
|
||||
|
||||
}
|
||||
if(pcevent.multi1==1 && pcevent.multi2==2) {
|
||||
//if(pcevent.multi1==1) {
|
||||
|
||||
bool TCC = sx3event.Time1 - cathodeT < 0;
|
||||
bool TCA = sx3event.Time1 - anodeT < 0;
|
||||
|
||||
plotter->Fill2D("pcz_vs_sx3pczguess_A1C2",600,-200,200,600,-200,200,pczguess,pcevent.pos.Z());
|
||||
double pcz_fix = pcfix_func.Eval(pcevent.pos.Z());
|
||||
|
||||
|
|
@ -969,49 +1026,14 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
plotter->Fill2D("phi_vs_stripnum",180,-180,180,48,0,48,pcevent.pos.Phi()*180./M_PI,sx3event.ch2);
|
||||
plotter->Fill2D("E_theta_AnodeSX3",400,-20,180,300,0,15,sx3theta*180/M_PI,sx3event.Energy1);
|
||||
}
|
||||
if(PCSX3TimeCut) {
|
||||
plotter->Fill1D("PCZ_sx3",800,-200,200,pcevent.pos.Z(),"hPCZSX3");
|
||||
plotter->Fill1D("PCZ",800,-200,200,pcevent.pos.Z(),"phicut");
|
||||
/*for(auto cc: cClusters)
|
||||
for(auto ac: aClusters) {
|
||||
plotter->Fill1D("PCZsx3_phicut_a"+std::to_string(ac.size())+"_c"+std::to_string(cc.size()),800,-200,200,pcevent.pos.Z(),"hPCZSX3");
|
||||
}*/
|
||||
}
|
||||
}//end PC-SX3 coincidence
|
||||
|
||||
/*for(size_t ii=0; ii<QQQ_Events.size(); ii++) {
|
||||
for(size_t jj=ii+1; jj<QQQ_Events.size(); jj++) {
|
||||
//if(TMath::Abs(QQQ_Events.at(ii).pos.Phi()*180/M_PI-QQQ_Events.at(jj).pos.Phi()*180/M_PI)>20) continue;
|
||||
if(QQQ_Events.at(ii).ch1 == QQQ_Events.at(jj).ch1) continue;
|
||||
if(QQQ_Events.at(ii).ch2 == QQQ_Events.at(jj).ch2) continue;
|
||||
if(QQQ_Events.at(ii).ch1 == QQQ_Events.at(jj).ch1-1) continue;
|
||||
if(QQQ_Events.at(ii).ch2 == QQQ_Events.at(jj).ch2-1) continue;
|
||||
if(QQQ_Events.at(ii).ch1 == QQQ_Events.at(jj).ch1+1) continue;
|
||||
if(QQQ_Events.at(ii).ch2 == QQQ_Events.at(jj).ch2+1) continue;
|
||||
for(auto qqqevent: QQQ_Events) {
|
||||
/* Events with QQQ present, but PC events don't have a reliable cathode signal, so we scan just the anode wire clusters */
|
||||
|
||||
double dt = QQQ_Events.at(ii).Time1-QQQ_Events.at(jj).Time1;
|
||||
plotter->Fill1D("dt_qqqi_qqqj",800,-2000,2000,dt);
|
||||
if(TMath::Abs(dt) > 150) continue;
|
||||
plotter->Fill1D("dt_qqqi_qqqj_coinc",800,-2000,2000,dt);
|
||||
double sum_e = QQQ_Events.at(ii).Energy1+QQQ_Events.at(jj).Energy1;
|
||||
plotter->Fill2D("sum_qqqE",400,0,30,400,0,30,QQQ_Events.at(ii).Energy1,sum_e);
|
||||
plotter->Fill2D("qqq_matrix",400,0,30,400,0,30,QQQ_Events.at(ii).Energy1,QQQ_Events.at(jj).Energy1);
|
||||
plotter->Fill2D("qqq_matrix",400,0,30,400,0,30,QQQ_Events.at(jj).Energy1,QQQ_Events.at(ii).Energy1);
|
||||
plotter->Fill2D("qqq_ch2_ch2",400,0,400,400,0,400,QQQ_Events.at(jj).ch2,QQQ_Events.at(ii).ch2);
|
||||
plotter->Fill2D("qqq_ch1_ch1",400,0,400,400,0,400,QQQ_Events.at(jj).ch1,QQQ_Events.at(ii).ch1);
|
||||
|
||||
if(sum_e > 6.50 && sum_e < 7.50) {
|
||||
plotter->Fill2D("qqq_ang1_ang2",180,-360,360,180,-360,360,QQQ_Events.at(jj).pos.Phi()*180/M_PI,QQQ_Events.at(ii).pos.Phi()*180/M_PI);
|
||||
//if(PC_Events.size()<2) continue;
|
||||
for(auto pcevent: PC_Events) {
|
||||
plotter->Fill2D("pcphi_vs_qqqphi_i_esumcut",180,-360,360,180,-360,360,pcevent.pos.Phi()*180/M_PI,QQQ_Events.at(ii).pos.Phi()*180/M_PI);
|
||||
plotter->Fill2D("pcphi_vs_qqqphi_j_esumcut",180,-360,360,180,-360,360,pcevent.pos.Phi()*180/M_PI,QQQ_Events.at(jj).pos.Phi()*180/M_PI);
|
||||
}
|
||||
}
|
||||
}
|
||||
}*/
|
||||
for(auto pcevent: PC_Events) {
|
||||
for(auto qqqevent: QQQ_Events) {
|
||||
|
||||
for(auto pcevent: PC_Events) {
|
||||
plotter->Fill1D("dt_pcA_qqqR",640,-2000,2000,qqqevent.Time1 - pcevent.Time1);
|
||||
plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE",640,-2000,2000,400,0,30,qqqevent.Time1-pcevent.Time1, qqqevent.Energy1);
|
||||
plotter->Fill1D("dt_pcC_qqqW",640,-2000,2000,qqqevent.Time2 - pcevent.Time2);
|
||||
|
|
@ -1086,6 +1108,16 @@ Bool_t MakeVertex::Process(Long64_t entry)
|
|||
double pcz_guess_int = z_to_crossover_rho(pcevent.pos.Z())/TMath::Tan((qqqevent.pos-TVector3(0,0,source_vertex)).Theta()) + source_vertex;
|
||||
//plotter->Fill2D("pczguess_vs_pc_int",180,0,200,150,0,200,pcz_guess_int,pcevent.pos.Z(),"phicut");
|
||||
plotter->Fill2D("pczguess_vs_pc_int",400,-200,200,600,-400,400,pcz_guess_int,pcevent.pos.Z(),"phicut");
|
||||
if(pcevent.multi1==1 && pcevent.multi2==1) {
|
||||
double pcz_fix_a1c1 = pcfix_func_a1c1.Eval(pcevent.pos.Z());
|
||||
TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix_a1c1);
|
||||
TVector3 v = x2f-x1;
|
||||
double t_minimum = -1.0*(x1.X()*v.X()+x1.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y());
|
||||
TVector3 r_rhoMin_fix = x1 + t_minimum*v;
|
||||
plotter->Fill1D("pczfix_A1C1_1d_qqq",600,-200,200,pcz_fix_a1c1);
|
||||
plotter->Fill2D("pczfix_vs_qqqpczguess_A1C1",600,-200,200,600,-200,200,pcz_guess_int,pcz_fix_a1c1);
|
||||
plotter->Fill2D("pczguess_vs_pc_int_A1C1",400,-200,200,600,-400,400,pcz_guess_int,pcevent.pos.Z(),"phicut");
|
||||
}
|
||||
if(pcevent.multi1==1 && pcevent.multi2==2) {
|
||||
double pcz_fix = pcfix_func.Eval(pcevent.pos.Z());
|
||||
TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix);
|
||||
|
|
@ -1567,8 +1599,7 @@ void MakeVertex::Terminate()
|
|||
}
|
||||
|
||||
|
||||
void protonAlphaHistograms(HistPlotter* plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events){
|
||||
|
||||
void protonAlphaHistograms(HistPlotter* plotter, std::vector<Event> QQQ_Events, std::vector<Event> SX3_Events, std::vector<Event> PC_Events) {
|
||||
//Sidetrack for a(p,p)
|
||||
std::string aplabel = "a(p,p)";
|
||||
Kinematics apkin_p(1.008664916,4.002603254,1.008664916,4.002603254,7.0);//m3 is proton
|
||||
|
|
@ -1586,7 +1617,7 @@ void protonAlphaHistograms(HistPlotter* plotter, std::vector<Event> QQQ_Events,
|
|||
|
||||
for(auto pcevent: PC_Events) {
|
||||
|
||||
double pcz_fix = pcfix_func.Eval(pcevent.pos.Z())-5.0;
|
||||
double pcz_fix = pcfix_func.Eval(pcevent.pos.Z())-15.0;
|
||||
TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix);
|
||||
TVector3 x1(qqqevent.pos);
|
||||
TVector3 v = x2f-x1;
|
||||
|
|
|
|||
845
MakeVertexSX3.C
845
MakeVertexSX3.C
|
|
@ -73,151 +73,36 @@ void MakeVertexSX3::Begin(TTree * /*tree*/)
|
|||
pw_contr.ConstructGeo();
|
||||
pwinstance.ConstructGeo();
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// 1. CRITICAL FIX: Initialize PC Arrays to Default (Raw)
|
||||
// ---------------------------------------------------------
|
||||
for (int i = 0; i < 48; i++)
|
||||
{
|
||||
pcSlope[i] = 1.0; // Default slope = 1 (preserves Raw energy)
|
||||
pcIntercept[i] = 0.0; // Default intercept = 0
|
||||
}
|
||||
|
||||
// Calculate Crossover Geometry ONCE
|
||||
TVector3 a, c, diff;
|
||||
double a2, ac, c2, adiff, cdiff, denom, alpha;
|
||||
|
||||
for (size_t i = 0; i < pwinstance.An.size(); i++)
|
||||
{
|
||||
a = pwinstance.An[i].first - pwinstance.An[i].second;
|
||||
|
||||
for (size_t j = 0; j < pwinstance.Ca.size(); j++)
|
||||
std::string dataset = "27Al";
|
||||
{
|
||||
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
|
||||
diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
|
||||
a2 = a.Dot(a);
|
||||
c2 = c.Dot(c);
|
||||
ac = a.Dot(c);
|
||||
adiff = a.Dot(diff);
|
||||
cdiff = c.Dot(diff);
|
||||
denom = a2 * c2 - ac * ac;
|
||||
alpha = (ac * cdiff - c2 * adiff) / denom;
|
||||
std::ifstream infile("sx3cal/"+dataset+"/backgains.dat");
|
||||
std::string temp;
|
||||
int backpos, frontpos, clkpos;
|
||||
if (infile.is_open())
|
||||
while(infile>>clkpos>>temp>>frontpos>>temp>>backpos>>sx3BackGain[clkpos][frontpos][backpos])
|
||||
;//std::cout << sx3BackGain[clkpos][frontpos][backpos] << std::endl;
|
||||
infile.close();
|
||||
|
||||
Crossover[i][j][0].x = pwinstance.An[i].first.X() + alpha * a.X();
|
||||
Crossover[i][j][0].y = pwinstance.An[i].first.Y() + alpha * a.Y();
|
||||
Crossover[i][j][0].z = pwinstance.An[i].first.Z() + alpha * a.Z();
|
||||
infile.open("sx3cal/"+dataset+"/frontgains.dat");
|
||||
if (infile.is_open())
|
||||
while(infile>>clkpos>>temp>>temp>>frontpos>>sx3FrontOffset[clkpos][frontpos]>>sx3FrontGain[clkpos][frontpos])
|
||||
;//std::cout << sx3FrontOffset[clkpos][frontpos] << " " << sx3FrontGain[clkpos][frontpos] << std::endl;
|
||||
infile.close();
|
||||
|
||||
if (Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190 || (i+j)%24 == 12)
|
||||
{
|
||||
Crossover[i][j][0].z = 9999999;
|
||||
}
|
||||
|
||||
Crossover[i][j][1].x = alpha;
|
||||
Crossover[i][j][1].y = 0;
|
||||
infile.open("sx3cal/"+dataset+"/rightgains.dat");
|
||||
if (infile.is_open())
|
||||
while(infile>>clkpos>>frontpos>>temp>>sx3RightGain[clkpos][frontpos]) {
|
||||
sx3RightGain[clkpos][frontpos]=TMath::Abs(sx3RightGain[clkpos][frontpos]);
|
||||
}
|
||||
infile.close();
|
||||
}
|
||||
}
|
||||
|
||||
// Load PC Calibrations
|
||||
std::ifstream inputFile("slope_intercept_results.txt");
|
||||
if (inputFile.is_open())
|
||||
{
|
||||
std::string line;
|
||||
int index;
|
||||
double slope, intercept;
|
||||
while (std::getline(inputFile, line))
|
||||
{
|
||||
std::stringstream ss(line);
|
||||
ss >> index >> slope >> intercept;
|
||||
if (index >= 0 && index <= 47)
|
||||
{
|
||||
pcSlope[index] = slope;
|
||||
pcIntercept[index] = intercept;
|
||||
}
|
||||
}
|
||||
inputFile.close();
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cerr << "Error opening slope_intercept.txt" << std::endl;
|
||||
}
|
||||
|
||||
// Load QQQ Cuts from file
|
||||
// {
|
||||
// std::string filename = "QQQ_PCCut.root";
|
||||
// TFile *cutFile = TFile::Open(filename.c_str(), "READ");
|
||||
// if (cutFile && !cutFile->IsZombie())
|
||||
// {
|
||||
// cutQQQ = (TCutg *)cutFile->Get("cutQQQPC");
|
||||
// if (cutQQQ)
|
||||
// {
|
||||
// std::cout << "Loaded QQQ PC cut from " << filename << std::endl;
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// std::cerr << "Error: cutQQQPC not found in " << filename << std::endl;
|
||||
// }
|
||||
// cutFile->Close();
|
||||
// }
|
||||
// }
|
||||
|
||||
// ... (Load QQQ Gains and Calibs - same as before) ...
|
||||
{
|
||||
std::string filename = "qqq_GainMatch.dat";
|
||||
std::ifstream infile(filename);
|
||||
if (infile.is_open())
|
||||
{
|
||||
int det, ring, wedge;
|
||||
double gainw, gainr;
|
||||
while (infile >> det >> wedge >> ring >> gainw >> gainr)
|
||||
{
|
||||
qqqGain[det][wedge][ring] = gainw;
|
||||
qqqGainValid[det][wedge][ring] = (gainw > 0);
|
||||
// std::cout << "QQQ Gain Loaded: Det " << det << " Ring " << ring << " Wedge " << wedge << " GainW " << gainw << " GainR " << gainr << std::endl;
|
||||
}
|
||||
infile.close();
|
||||
}
|
||||
}
|
||||
{
|
||||
std::string filename = "qqq_Calib.dat";
|
||||
std::ifstream infile(filename);
|
||||
if (infile.is_open())
|
||||
{
|
||||
int det, ring, wedge;
|
||||
double slope;
|
||||
while (infile >> det >> wedge >> ring >> slope)
|
||||
{
|
||||
qqqCalib[det][wedge][ring] = slope;
|
||||
qqqCalibValid[det][wedge][ring] = (slope > 0);
|
||||
// std::cout << "QQQ Calib Loaded: Det " << det << " Ring " << ring << " Wedge " << wedge << " Slope " << slope << std::endl;
|
||||
}
|
||||
infile.close();
|
||||
}
|
||||
}
|
||||
|
||||
{
|
||||
std::ifstream infile("sx3cal/backgains.dat");
|
||||
std::string temp;
|
||||
int backpos, frontpos, clkpos;
|
||||
std::cout << "foo" << std::endl;
|
||||
if (infile.is_open())
|
||||
while(infile>>clkpos>>temp>>frontpos>>temp>>backpos>>sx3BackGain[clkpos][frontpos][backpos])
|
||||
std::cout << sx3BackGain[clkpos][frontpos][backpos] << std::endl;
|
||||
infile.close();
|
||||
|
||||
infile.open("sx3cal/frontgains.dat");
|
||||
if (infile.is_open())
|
||||
while(infile>>clkpos>>temp>>temp>>frontpos>>sx3FrontOffset[clkpos][frontpos]>>sx3FrontGain[clkpos][frontpos])
|
||||
std::cout << sx3FrontOffset[clkpos][frontpos] << " " << sx3FrontGain[clkpos][frontpos] << std::endl;
|
||||
infile.close();
|
||||
|
||||
infile.open("sx3cal/rightgains.dat");
|
||||
if (infile.is_open())
|
||||
while(infile>>clkpos>>frontpos>>temp>>sx3RightGain[clkpos][frontpos]) {
|
||||
sx3RightGain[clkpos][frontpos]=TMath::Abs(sx3RightGain[clkpos][frontpos]);
|
||||
}
|
||||
infile.close();
|
||||
|
||||
}
|
||||
std::cout << "aaa" << std::endl;
|
||||
}
|
||||
|
||||
Bool_t MakeVertexSX3::Process(Long64_t entry)
|
||||
|
|
@ -257,7 +142,7 @@ Bool_t MakeVertexSX3::Process(Long64_t entry)
|
|||
if(sx3.ch[i]>=8) {
|
||||
int sx3ch=sx3.ch[i]-8;
|
||||
sx3ch=(sx3ch+3)%4;
|
||||
if(sx3ch==0 || sx3ch==3) continue;
|
||||
//if(sx3ch==0 || sx3ch==3) continue;
|
||||
float value=sx3.e[i];
|
||||
int gch = sx3.id[i]*4+(sx3.ch[i]-8);
|
||||
Fsx3.at(id).fillevent("BACK",sx3ch,value);
|
||||
|
|
@ -293,666 +178,44 @@ Bool_t MakeVertexSX3::Process(Long64_t entry)
|
|||
Event sx3ev(TVector3(0,0,z),backE,-1,det.ts,-1,det.stripB+4*id,det.stripF+4*id);
|
||||
sx3Events.push_back(sx3ev);
|
||||
}
|
||||
}
|
||||
}
|
||||
//return kTRUE;
|
||||
// QQQ Processing
|
||||
|
||||
int qqqCount = 0;
|
||||
int qqqAdjCh = 0;
|
||||
// REMOVE WHEN RERUNNING USING THE NEW CALIBRATION FILE
|
||||
for (int i = 0; i < qqq.multi; i++)
|
||||
{
|
||||
//if ((qqq.id[i] == 3 || qqq.id[i] == 1) && qqq.ch[i] < 16)
|
||||
if (qqq.id[i] == 1 && qqq.ch[i] < 16) //for run 12, 26Al
|
||||
{
|
||||
qqq.ch[i] = 16 - qqq.ch[i];
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < qqq.multi; i++)
|
||||
{
|
||||
if (qqq.id[i] == 0 && qqq.ch[i] >= 16)
|
||||
{
|
||||
qqq.ch[i] = 31 - qqq.ch[i] + 16;
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<Event> QQQ_Events, PC_Events;
|
||||
std::vector<Event> QQQ_Events_Raw, PC_Events_Raw;
|
||||
std::vector<Event> QQQ_Events2; //clustering done
|
||||
|
||||
std::unordered_map<int,std::tuple<int,int,double,double>> qvecr[4], qvecw[4];
|
||||
if(qqq.multi>1) {
|
||||
//if(qqq.multi>=3) std::cout << "-----" << std::endl;
|
||||
for(int i=0; i<qqq.multi; i++) {
|
||||
//if(qqq.multi>=3) std::cout << std::setprecision(16) << "qqq"<< qqq.id[i] << " " << std::string(qqq.ch[i]/16?"ring":"wedge") << qqq.ch[i]%16 << " " << qqq.e[i] << " " << qqq.t[i] - qqq.t[0] << std::endl;
|
||||
if(qqq.ch[i]/16) {
|
||||
if(qvecr[qqq.id[i]].find(qqq.ch[i])!=qvecr[qqq.id[i]].end()) std::cout << "mayday!" << std::endl;
|
||||
qvecr[qqq.id[i]][qqq.ch[i]] = std::tuple(qqq.id[i],qqq.ch[i],qqq.e[i],qqq.t[i]);
|
||||
} else {
|
||||
if(qvecw[qqq.id[i]].find(qqq.ch[i])!=qvecw[qqq.id[i]].end()) std::cout << "mayday!" << std::endl;
|
||||
qvecw[qqq.id[i]][qqq.ch[i]] = std::tuple(qqq.id[i],qqq.ch[i],qqq.e[i],qqq.t[i]);
|
||||
if(!det.valid && det.valid_back_chans.size()) {
|
||||
for(auto bch : det.valid_back_chans) {
|
||||
plotter->Fill1D("good_backs_bad_fronts_id"+std::to_string(id)+"_b"+std::to_string(bch),800,0,4096,det.back[bch].at(0),"good_backs_no_front");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
bool PCQQQTimeCut = false;
|
||||
for (int i = 0; i < qqq.multi; i++) {
|
||||
plotter->Fill2D("QQQ_Index_Vs_Energy", 16 * 8, 0, 16 * 8, 2000, 0, 16000, qqq.index[i], qqq.e[i], "hRawQQQ");
|
||||
|
||||
for (int j = 0; j < qqq.multi; j++) {
|
||||
if (j == i)
|
||||
continue;
|
||||
plotter->Fill2D("QQQ_Coincidence_Matrix", 16 * 8, 0, 16 * 8, 16 * 8, 0, 16 * 8, qqq.index[i], qqq.index[j], "hRawQQQ");
|
||||
}
|
||||
|
||||
for (int k = 0; k < pc.multi; k++) {
|
||||
if (pc.index[k] < 24 && pc.e[k] > 50) {
|
||||
plotter->Fill2D("QQQ_Vs_Anode_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k], "hRawQQQ");
|
||||
plotter->Fill2D("QQQ_Vs_PC_Index", 16 * 8, 0, 16 * 8, 24, 0, 24, qqq.index[i], pc.index[k], "hRawQQQ");
|
||||
}
|
||||
else if (pc.index[k] >= 24 && pc.e[k] > 50) {
|
||||
plotter->Fill2D("QQQ_Vs_Cathode_Energy", 400, 0, 4000, 1000, 0, 16000, qqq.e[i], pc.e[k], "hRawQQQ");
|
||||
}
|
||||
}
|
||||
|
||||
for (int j = i + 1; j < qqq.multi; j++) {
|
||||
if (qqq.id[i] == qqq.id[j]) {
|
||||
qqqCount++;
|
||||
|
||||
int chWedge = -1;
|
||||
int chRing = -1;
|
||||
double eWedge = 0.0;
|
||||
double eWedgeMeV = 0.0;
|
||||
double eRing = 0.0;
|
||||
double eRingMeV = 0.0;
|
||||
double tRing = 0.0;
|
||||
double tWedge = 0.0;
|
||||
|
||||
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]) {
|
||||
chWedge = qqq.ch[i];
|
||||
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
|
||||
chRing = qqq.ch[j] - 16;
|
||||
eRing = qqq.e[j];
|
||||
tRing = static_cast<double>(qqq.t[j]);
|
||||
tWedge = static_cast<double>(qqq.t[i]);
|
||||
}
|
||||
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]) {
|
||||
chWedge = qqq.ch[j];
|
||||
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
|
||||
chRing = qqq.ch[i] - 16;
|
||||
eRing = qqq.e[i];
|
||||
tRing = static_cast<double>(qqq.t[i]);
|
||||
tWedge = static_cast<double>(qqq.t[j]);
|
||||
}
|
||||
else
|
||||
continue;
|
||||
|
||||
plotter->Fill1D("Wedgetime_Vs_Ringtime", 100, -1000, 1000, tWedge - tRing, "hTiming");
|
||||
plotter->Fill2D("RingE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chRing + qqq.id[i] * 16, eRing, "hRawQQQ");
|
||||
plotter->Fill2D("WedgeE_vs_Index", 16 * 4, 0, 16 * 4, 1000, 0, 16000, chWedge + qqq.id[i] * 16, eWedge, "hRawQQQ");
|
||||
|
||||
if (qqqCalibValid[qqq.id[i]][chWedge][chRing]) {
|
||||
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000;
|
||||
eRingMeV = eRing * qqqCalib[qqq.id[i]][chWedge][chRing] / 1000;
|
||||
|
||||
if(eRingMeV/eWedgeMeV > 3.0 || eRingMeV/eWedgeMeV<1.0/3.0) continue;
|
||||
//if(eRingMeV<4.0 || eWedgeMeV<4.0) continue;
|
||||
|
||||
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
|
||||
double rho = 50. + (50. / 16.) * (chRing + 0.5); //"?"
|
||||
//z used to be 75+30+23=128
|
||||
//we found a 12mm shift towards the vertex later --> 116
|
||||
Event qqqevent(TVector3(rho*TMath::Cos(theta),rho*TMath::Sin(theta),116), eRingMeV, eWedgeMeV, tRing, tWedge,chRing+qqq.id[i]*16, chWedge+qqq.id[i]*16);
|
||||
Event qqqeventr(TVector3(rho*TMath::Cos(theta),rho*TMath::Sin(theta),116), eRing, eWedge, tRing, tWedge,chRing+qqq.id[i]*16, chWedge+qqq.id[i]*16);
|
||||
QQQ_Events.push_back(qqqevent);
|
||||
QQQ_Events_Raw.push_back(qqqeventr);
|
||||
plotter->Fill2D("QQQCartesianPlot", 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
|
||||
plotter->Fill2D("QQQCartesianPlot" + std::to_string(qqq.id[i]), 200, -100, 100, 200, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hCalQQQ");
|
||||
if (PCQQQTimeCut) {
|
||||
plotter->Fill2D("PC_XY_Projection_QQQ_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
|
||||
}
|
||||
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, rho * TMath::Cos(theta), rho * TMath::Sin(theta), "hPCQQQ");
|
||||
}
|
||||
else
|
||||
continue;
|
||||
|
||||
plotter->Fill2D("WedgeE_Vs_RingECal", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
|
||||
plotter->Fill2D("WedgeE_Vs_RingECal_selected", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ");
|
||||
|
||||
for (int k = 0; k < pc.multi; k++)
|
||||
{
|
||||
plotter->Fill2D("RingCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
|
||||
plotter->Fill2D("WedgeCh_vs_Anode_Index", 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
|
||||
plotter->Fill2D("WedgeCh_vs_Anode_Index" + std::to_string(qqq.id[i]), 16 * 4, 0, 16 * 4, 24, 0, 24, chWedge + qqq.id[i] * 16, pc.index[k]);
|
||||
plotter->Fill2D("RingCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chRing + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
|
||||
plotter->Fill2D("WedgeCh_vs_Cathode_Index", 16 * 4, 0, 16 * 4, 24, 24, 48, chWedge + qqq.id[i] * 16, pc.index[k], "hRawQQQ");
|
||||
|
||||
if (pc.index[k] < 24 && pc.e[k] > 50)
|
||||
{
|
||||
plotter->Fill2D("Timing_Difference_QQQ_PC", 500, -2000, 2000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming");
|
||||
plotter->Fill2D("DelT_Vs_QQQRingECal", 500, -2000, 2000, 1000, 0, 10, tRing - static_cast<double>(pc.t[k]), eRingMeV, "hTiming");
|
||||
plotter->Fill2D("CalibratedQQQEvsPCE_R", 1000, 0, 10, 2000, 0, 30000, eRingMeV, pc.e[k], "hPCQQQ");
|
||||
plotter->Fill2D("CalibratedQQQEvsPCE_W", 1000, 0, 10, 2000, 0, 30000, eWedgeMeV, pc.e[k], "hPCQQQ");
|
||||
if (tRing - static_cast<double>(pc.t[k]) < -150) // proton tests, 27Al
|
||||
//if (tRing - static_cast<double>(pc.t[k]) < -150 && tRing - static_cast<double>(pc.t[k]) > -450) // 27Al
|
||||
//if (tRing - static_cast<double>(pc.t[k]) < -70 && tRing - static_cast<double>(pc.t[k]) > -150) // 17F
|
||||
{
|
||||
PCQQQTimeCut = true;
|
||||
}
|
||||
}
|
||||
|
||||
if (pc.index[k] >= 24 && pc.e[k] > 50) {
|
||||
plotter->Fill2D("Timing_Difference_QQQ_PC_Cathode", 500, -2000, 2000, 16, 0, 16, tRing - static_cast<double>(pc.t[k]), chRing, "hTiming");
|
||||
}
|
||||
} //end of pc k loop
|
||||
|
||||
if (!HitNonZero) {
|
||||
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
|
||||
double rho = 50. + (50. / 16.) * (chRing + 0.5); //"?"
|
||||
double x = rho * TMath::Cos(theta);
|
||||
double y = rho * TMath::Sin(theta);
|
||||
hitPos.SetXYZ(x, y, (23 + 75 + 30));
|
||||
qqqenergy = eRingMeV;
|
||||
qqqtimestamp = tRing;
|
||||
HitNonZero = true;
|
||||
}
|
||||
} // if j==i
|
||||
} //j loop end
|
||||
} //i loop end
|
||||
|
||||
plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hRawQQQ");
|
||||
|
||||
/*if(QQQ_Events.size()>=1) {
|
||||
std::cout<< " ---->" << std::endl;
|
||||
for(auto qe: QQQ_Events) {
|
||||
std::cout << qe.ch1/16 << " " <<qe.ch2/16 << " " << qe.ch1%16 << " "<< qe.ch2%16 << " " << qe.Energy1 << " " << qe.Energy2 << " " << std::endl;
|
||||
}
|
||||
}*/
|
||||
|
||||
|
||||
typedef std::unordered_map<int,std::tuple<int,double,double>> WireEvent; //this stores nearest neighbour wire events, or a 'cluster'
|
||||
WireEvent aWireEvents, cWireEvents; //naming for book keeping
|
||||
aWireEvents.clear();
|
||||
aWireEvents.reserve(24);
|
||||
|
||||
// PC Gain Matching and Filling
|
||||
double anodeT = -99999;
|
||||
double cathodeT = 99999;
|
||||
int anodeIndex = -1;
|
||||
int cathodeIndex = -1;
|
||||
for (int i = 0; i < pc.multi; i++)
|
||||
{
|
||||
if (pc.e[i] > 50)
|
||||
{
|
||||
plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], static_cast<double>(pc.e[i]), "hRawPC");
|
||||
} else {
|
||||
continue;
|
||||
}
|
||||
|
||||
if (pc.index[i] < 48)
|
||||
{
|
||||
pc.e[i] = pcSlope[pc.index[i]] * pc.e[i] + pcIntercept[pc.index[i]];
|
||||
plotter->Fill2D("PC_Index_VS_GainMatched_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], pc.e[i], "hGMPC");
|
||||
}
|
||||
|
||||
if (pc.index[i] < 24)
|
||||
{
|
||||
anodeT = static_cast<double>(pc.t[i]);
|
||||
anodeIndex = pc.index[i];
|
||||
aWireEvents[pc.index[i]] = std::tuple(pc.index[i],pc.e[i],static_cast<double>(pc.t[i]));
|
||||
}
|
||||
else
|
||||
{
|
||||
cathodeT = static_cast<double>(pc.t[i]);
|
||||
cathodeIndex = pc.index[i] - 24;
|
||||
cWireEvents[pc.index[i]-24] = std::tuple(pc.index[i]-24,pc.e[i],static_cast<double>(pc.t[i]));
|
||||
}
|
||||
|
||||
if (anodeT != -99999 && cathodeT != 99999)
|
||||
{
|
||||
for (int j = 0; j < qqq.multi; j++)
|
||||
{
|
||||
plotter->Fill1D("PC_Time_qqq", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
|
||||
plotter->Fill2D("PC_Time_Vs_QQQ_ch", 200, -2000, 2000, 16 * 8, 0, 16 * 8, anodeT - cathodeT, qqq.ch[j], "hTiming");
|
||||
plotter->Fill2D("PC_Time_vs_AIndex", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, anodeIndex, "hTiming");
|
||||
plotter->Fill2D("PC_Time_vs_CIndex", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, cathodeIndex, "hTiming");
|
||||
// plotter->Fill1D("PC_Time_A" + std::to_string(anodeIndex) + "_C" + std::to_string(cathodeIndex), 200, -1000, 1000, anodeT - cathodeT, "TimingPC");
|
||||
}
|
||||
|
||||
for (int j = 0; j < sx3.multi; j++)
|
||||
{
|
||||
plotter->Fill1D("PC_Time_sx3", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
|
||||
}
|
||||
|
||||
plotter->Fill1D("PC_Time", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
|
||||
}
|
||||
|
||||
for (int j = i + 1; j < pc.multi; j++)
|
||||
{
|
||||
plotter->Fill2D("PC_Coincidence_Matrix", 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
|
||||
plotter->Fill2D("PC_Coincidence_Matrix_anodeMinusCathode_lt_-200_" + std::to_string(anodeT - cathodeT < -200), 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
|
||||
plotter->Fill2D("Anode_V_Anode", 24, 0, 24, 24, 0, 24, pc.index[i], pc.index[j], "hGMPC");
|
||||
}
|
||||
}
|
||||
|
||||
anodeHits.clear();
|
||||
cathodeHits.clear();
|
||||
corrcatMax.clear();
|
||||
|
||||
int aID = 0;
|
||||
int cID = 0;
|
||||
double aE = 0;
|
||||
double cE = 0;
|
||||
double aESum = 0;
|
||||
double cESum = 0;
|
||||
double aEMax = 0;
|
||||
int aIDMax = 0;
|
||||
|
||||
|
||||
|
||||
for (int i = 0; i < pc.multi; i++) {
|
||||
// if (pc.e[i] > 100)
|
||||
{
|
||||
if (pc.index[i] < 24) {
|
||||
anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
|
||||
}
|
||||
else if (pc.index[i] >= 24) {
|
||||
cathodeHits.push_back(std::pair<int, double>(pc.index[i] - 24, pc.e[i]));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::sort(anodeHits.begin(),anodeHits.end(),[](std::pair<int,double> a, std::pair<int,double> b){ return a.first < b.first;});
|
||||
std::sort(cathodeHits.begin(),cathodeHits.end(),[](std::pair<int,double> a, std::pair<int,double> b){ return a.first < b.first;});
|
||||
|
||||
//clusters = collection of (collection of wires) where each wire is (index, energy, timestamp)
|
||||
std::vector<std::vector<std::tuple<int,double,double>>> aClusters = pwinstance.Make_Clusters(aWireEvents);
|
||||
std::vector<std::vector<std::tuple<int,double,double>>> cClusters = pwinstance.Make_Clusters(cWireEvents);
|
||||
|
||||
std::vector<std::pair<double,double>> sumE_AC;
|
||||
for(auto aCluster: aClusters) {
|
||||
for(auto cCluster: cClusters) {
|
||||
if(aCluster.size()<=1 && cCluster.size()<=1) continue;
|
||||
auto [crossover,alpha,apSumE,cpSumE,apMaxE,cpMaxE,apTSMaxE,cpTSMaxE] = pwinstance.FindCrossoverProperties(aCluster, cCluster);
|
||||
if(alpha!=9999999 && apSumE!=-1) {
|
||||
//Event PCEvent(crossover,apMaxE,cpMaxE,apTSMaxE,cpTSMaxE);
|
||||
//Event PCEvent(crossover,apSumE,cpSumE,apTSMaxE,cpTSMaxE);
|
||||
Event PCEvent(crossover,apSumE,cpMaxE,apTSMaxE,cpTSMaxE); //run12 shows cathode-max and anode-sum provide best dE signals.
|
||||
//std::cout << apSumE << " " << crossover.Perp() << " " << apMaxE << " " << apTSMaxE << std::endl;
|
||||
PC_Events.push_back(PCEvent);
|
||||
sumE_AC.push_back(std::pair(apSumE,cpSumE));
|
||||
}
|
||||
}
|
||||
}
|
||||
if(QQQ_Events.size() && PC_Events.size())
|
||||
plotter->Fill2D("PCEv_vs_QQQEv",20,0,20,20,0,20,QQQ_Events.size(),PC_Events.size());
|
||||
|
||||
for(auto pcevent:PC_Events) {
|
||||
for(auto sx3event:sx3Events) {
|
||||
plotter->Fill1D("dt_pcA_sx3B"+std::to_string(sx3event.ch2),640,-2000,2000,sx3event.Time1 - pcevent.Time1);
|
||||
plotter->Fill1D("dt_pcC_sx3B"+std::to_string(sx3event.ch2),640,-2000,2000,sx3event.Time1 - pcevent.Time2);
|
||||
plotter->Fill2D("dE_E_Anodesx3B",400,0,10,800,0,40000,sx3event.Energy1*0.001,pcevent.Energy1);
|
||||
|
||||
plotter->Fill2D("dE_E_Cathodesx3B",400,0,10,800,0,10000,sx3event.Energy1*0.001,pcevent.Energy2);
|
||||
double sx3z = sx3event.pos.Z()+(75.0/2.0)+23.0-90.0; //w.r.t target origin at 90 for run12
|
||||
double sx3rho = 88.0;
|
||||
double sx3theta = TMath::ATan2(sx3rho,sx3z);
|
||||
double pczguess = 40.0/TMath::Tan(sx3theta) + 90.0;
|
||||
plotter->Fill2D("pcz_vs_sx3pczguess",300,0,200,150,0,200,pczguess,pcevent.pos.Z());
|
||||
plotter->Fill2D("pcz_vs_sx3pczguess"+std::to_string(sx3event.ch2),300,0,200,150,0,200,pczguess,pcevent.pos.Z());
|
||||
plotter->Fill2D("pcz_vs_sx3z",300,0,200,150,0,200,sx3z+90,pcevent.pos.Z());
|
||||
}
|
||||
}
|
||||
|
||||
for(auto pcevent: PC_Events) {
|
||||
for(auto qqqevent: QQQ_Events) {
|
||||
plotter->Fill1D("dt_pcA_qqqR",640,-2000,2000,qqqevent.Time1 - pcevent.Time1);
|
||||
plotter->Fill1D("dt_pcC_qqqW",640,-2000,2000,qqqevent.Time2 - pcevent.Time2);
|
||||
plotter->Fill2D("dE_E_AnodeQQQR",400,0,10,800,0,40000,qqqevent.Energy1,pcevent.Energy1);
|
||||
plotter->Fill2D("dE_E_CathodeQQQR",400,0,10,800,0,10000,qqqevent.Energy2,pcevent.Energy2);
|
||||
double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0,0,90)).Theta())/TMath::Sin((TVector3(51.5,0,128.) - TVector3(0,0,90)).Theta());
|
||||
plotter->Fill2D("dE2_E_AnodeQQQR",400,0,10,800,0,40000,qqqevent.Energy1,pcevent.Energy1*sinTheta);
|
||||
plotter->Fill2D("dE2_E_CathodeQQQR",400,0,10,800,0,10000,qqqevent.Energy2,pcevent.Energy2*sinTheta);
|
||||
|
||||
if(qqqevent.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/4. && qqqevent.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/4.) {
|
||||
plotter->Fill1D("PCZ",800,-200,200,pcevent.pos.Z(),"phicut");
|
||||
double pcz_guess = 40.0/TMath::Tan((qqqevent.pos-TVector3(0,0,90)).Theta()) + 90; //this is ideally kept to be all QQQ+userinput for calibration of pcz
|
||||
plotter->Fill2D("pczguess_vs_pc",300,0,200,150,0,200,pcz_guess,pcevent.pos.Z(),"phicut");
|
||||
plotter->Fill2D("pczguess_vs_pc_phi="+std::to_string(qqqevent.pos.Phi()*180./M_PI),300,0,200,150,0,200,pcz_guess,pcevent.pos.Z(),"phicut");
|
||||
//plotter->Fill1D("PCZ",800,-200,200,pcevent.pos.Z(),"phicut");
|
||||
}
|
||||
}
|
||||
}
|
||||
//HALFTIME! Can stop here in future versions
|
||||
//return kTRUE;
|
||||
|
||||
if (anodeHits.size() >= 1 && cathodeHits.size() >= 1)
|
||||
{
|
||||
// 2. CRITICAL FIX: Define reference vector 'a'
|
||||
// In Analyzer.cxx, 'a' was left over from the loop. We use the first anode wire as reference here.
|
||||
// (Assuming pwinstance.An is populated and wires are generally parallel).
|
||||
TVector3 refAnode = pwinstance.An[0].first - pwinstance.An[0].second;
|
||||
|
||||
{
|
||||
for (const auto &anode : anodeHits)
|
||||
{
|
||||
aID = anode.first;
|
||||
aE = anode.second;
|
||||
aESum += aE;
|
||||
if (aE > aEMax)
|
||||
{
|
||||
aEMax = aE;
|
||||
aIDMax = aID;
|
||||
}
|
||||
}
|
||||
|
||||
for (const auto &cathode : cathodeHits)
|
||||
{
|
||||
cID = cathode.first;
|
||||
cE = cathode.second;
|
||||
plotter->Fill2D("AnodeMax_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aIDMax, cID, "hRawPC");
|
||||
plotter->Fill2D("Anode_Vs_Cathode_Coincidence_Matrix", 24, 0, 24, 24, 0, 24, aID, cID, "hRawPC");
|
||||
plotter->Fill2D("Anode_vs_CathodeE", 2000, 0, 30000, 2000, 0, 30000, aE, cE, "hGMPC");
|
||||
plotter->Fill2D("CathodeMult_V_CathodeE", 6, 0, 6, 2000, 0, 30000, cathodeHits.size(), cE, "hGMPC");
|
||||
for (int j = -4; j < 3; j++)
|
||||
{
|
||||
if ((aIDMax + 24 + j) % 24 == 23 - cID)
|
||||
{
|
||||
corrcatMax.push_back(std::pair<int, double>(cID, cE));
|
||||
cESum += cE;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
TVector3 anodeIntersection,vector_closest_to_z;
|
||||
anodeIntersection.Clear();
|
||||
vector_closest_to_z.Clear();
|
||||
if (corrcatMax.size() > 0)
|
||||
{
|
||||
double x = 0, y = 0, z = 0;
|
||||
for (const auto &corr : corrcatMax)
|
||||
{
|
||||
if (Crossover[aIDMax][corr.first][0].z > 9000000)
|
||||
continue;
|
||||
if (cESum > 0)
|
||||
{
|
||||
x += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].x;
|
||||
y += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].y;
|
||||
z += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].z;
|
||||
}
|
||||
}
|
||||
if (x == 0 && y == 0 && z == 0)
|
||||
;
|
||||
// to ignore events with no valid crossover points
|
||||
else
|
||||
anodeIntersection = TVector3(x, y, z);
|
||||
// << "Anode Intersection: " << anodeIntersection.X() << ", " << anodeIntersection.Y() << ", " << anodeIntersection.Z() << std::endl;
|
||||
}
|
||||
bool PCQQQPhiCut = false;
|
||||
// flip the algorithm for cathode 1 multi anode events
|
||||
if ((hitPos.Phi() > (anodeIntersection.Phi() - TMath::PiOver4())) && (hitPos.Phi() < (anodeIntersection.Phi() + TMath::PiOver4()))) {
|
||||
PCQQQPhiCut = true;
|
||||
}
|
||||
|
||||
for (double Tz = 60; Tz <= 100; Tz += 1.0)
|
||||
{
|
||||
TVector3 TargetPos(0, 0, Tz);
|
||||
if(PCQQQPhiCut && anodeIntersection.Perp()>0 && anodeIntersection.Z()!=0 && cathodeHits.size()>=2) {
|
||||
plotter->Fill2D("Inttheta_vs_QQQtheta_TC" + std::to_string(PCQQQTimeCut) + "_TZ" + std::to_string(Tz), 400, 0, 180, 90, 0, 90, (anodeIntersection - TargetPos).Theta() * 180. / TMath::Pi(), (hitPos - TargetPos).Theta() * 180. / TMath::Pi(), "TPosVariation");
|
||||
//plotter->Fill2D("R_ratio_to_Z_ratio" + std::to_string(PCQQQTimeCut) + "_TZ" + std::to_string(Tz), 100, -2, 2, 100, -2, 2, (anodeIntersection - TargetPos).Z()/(hitPos-TargetPos).Z(), ((anodeIntersection - TargetPos).Perp()+2.5)/(hitPos-TargetPos).Perp(), "TPosVariation");
|
||||
}
|
||||
}
|
||||
|
||||
if (anodeIntersection.Z() != 0 && anodeIntersection.Perp()>0 && HitNonZero)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_Projection", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
plotter->Fill2D("Z_Proj_VsDelTime", 600, -300, 300, 200, -2000, 2000, anodeIntersection.Z(), anodeT - cathodeT, "hPCzQQQ");
|
||||
plotter->Fill2D("IntPhi_vs_QQQphi", 100, -200, 200, 80, -200, 200, anodeIntersection.Phi() * 180. / TMath::Pi(), hitPos.Phi() * 180. / TMath::Pi(), "hPCQQQ");
|
||||
//plotter->Fill2D("Inttheta_vs_QQQtheta", 90, 0, 180, 20, 0, 45, anodeIntersection.Theta() * 180. / TMath::Pi(), hitPos.Theta() * 180. / TMath::Pi(), "hPCQQQ");
|
||||
//plotter->Fill2D("Inttheta_vs_QQQtheta_TC" + std::to_string(PCQQQTimeCut)+ "_PC"+std::to_string(PCQQQPhiCut), 90, 0, 180, 20, 0, 45, anodeIntersection.Theta() * 180. / TMath::Pi(), hitPos.Theta() * 180. / TMath::Pi(), "hPCQQQ");
|
||||
plotter->Fill2D("IntPhi_vs_QQQphi_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 100, -200, 200, 80, -200, 200, anodeIntersection.Phi() * 180. / TMath::Pi(), hitPos.Phi() * 180. / TMath::Pi(), "hPCQQQ");
|
||||
}
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() >= 2)
|
||||
plotter->Fill1D("PC_Z_Projection_TC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 1)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_1C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_1C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hPCzQQQ");
|
||||
}
|
||||
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_2C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_2C", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
|
||||
}
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() > 2)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_nC", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
plotter->Fill2D("IntersectionPhi_vs_AnodeZ_nC", 400, -200, 200, 600, -300, 300, anodeIntersection.Phi() * 180. / TMath::Pi(), anodeIntersection.Z(), "hGMPC");
|
||||
}
|
||||
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
|
||||
plotter->Fill2D("AHits_vs_CHits", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
|
||||
// make another plot with nearest neighbour constraint
|
||||
bool hasNeighbourAnodes = false;
|
||||
bool hasNeighbourCathodes = false;
|
||||
|
||||
// 1. Check Anodes for neighbours (including wrap-around 0-23)
|
||||
for (size_t i = 0; i < anodeHits.size(); i++)
|
||||
{
|
||||
for (size_t j = i + 1; j < anodeHits.size(); j++)
|
||||
{
|
||||
int diff = std::abs(anodeHits[i].first - anodeHits[j].first);
|
||||
if (diff == 1 || diff == 23)
|
||||
{ // 23 handles the cylindrical wrap
|
||||
hasNeighbourAnodes = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (hasNeighbourAnodes)
|
||||
break;
|
||||
}
|
||||
|
||||
// 2. Check Cathodes for neighbours (including wrap-around 0-23)
|
||||
for (size_t i = 0; i < cathodeHits.size(); i++)
|
||||
{
|
||||
for (size_t j = i + 1; j < cathodeHits.size(); j++)
|
||||
{
|
||||
int diff = std::abs(cathodeHits[i].first - cathodeHits[j].first);
|
||||
if (diff == 1 || diff == 23)
|
||||
{
|
||||
hasNeighbourCathodes = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (hasNeighbourCathodes)
|
||||
break;
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------
|
||||
// FILL PLOTS
|
||||
// ---------------------------------------------------------
|
||||
if (anodeHits.size() > 0 && cathodeHits.size() > 0)
|
||||
{
|
||||
plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
|
||||
// Constraint Plot: Only fill if BOTH planes have adjacent hits
|
||||
// This effectively removes events with only isolated single-wire hits (noise)
|
||||
if (hasNeighbourAnodes && hasNeighbourCathodes)
|
||||
{
|
||||
plotter->Fill2D("AHits_vs_CHits_NN", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
|
||||
}
|
||||
}
|
||||
|
||||
if (HitNonZero && anodeIntersection.Z() != 0)
|
||||
{
|
||||
pw_contr.CalTrack2(hitPos, anodeIntersection);
|
||||
plotter->Fill1D("VertexRecon", 600, -1300, 1300, pw_contr.GetZ0());
|
||||
plotter->Fill1D("VertexRecon_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, pw_contr.GetZ0());
|
||||
|
||||
if (cathodeHits.size() == 2)
|
||||
plotter->Fill1D("VertexRecon_2c_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, pw_contr.GetZ0());
|
||||
|
||||
TVector3 x2(anodeIntersection), x1(hitPos);
|
||||
|
||||
TVector3 v = x2-x1;
|
||||
double t_minimum = -1.0*(x1.X()*v.X()+x1.Y()*v.Y())/(v.X()*v.X()+v.Y()*v.Y());
|
||||
vector_closest_to_z = x1 + t_minimum*v;
|
||||
|
||||
plotter->Fill1D("VertexRecon_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
|
||||
if(vector_closest_to_z.Perp() < 20) {
|
||||
plotter->Fill1D("VertexRecon_RadialCut_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
|
||||
}
|
||||
|
||||
plotter->Fill2D("VertexRecon_XY_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 100, -100, 100, 100,-100,100, vector_closest_to_z.X(), vector_closest_to_z.Y() ,"customVertex");
|
||||
if(cathodeHits.size()==2) {
|
||||
plotter->Fill1D("VertexRecon2C_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
|
||||
if(vector_closest_to_z.Perp() < 20) {
|
||||
plotter->Fill1D("VertexRecon2C_RadialCut_Z_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, vector_closest_to_z.Z() ,"customVertex");
|
||||
}
|
||||
plotter->Fill2D("VertexRecon2C_XY_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 100, -100, 100, 100,-100,100, vector_closest_to_z.X(), vector_closest_to_z.Y() ,"customVertex");
|
||||
plotter->Fill2D("VertexRecon2C_RhoZ_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 100, -100, 100, 600,-1300,1300, vector_closest_to_z.Perp(), vector_closest_to_z.Z() ,"customVertex");
|
||||
plotter->Fill2D("VertexRecon2C_Z_vs_QQQE_TC"+std::to_string(PCQQQTimeCut)+"_PhiC"+std::to_string(PCQQQPhiCut), 600, -1300, 1300, 800,0,20000, vector_closest_to_z.Z(), qqqenergy ,"customVertex");
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
for (int i = 0; i < qqq.multi; i++)
|
||||
{
|
||||
if(anodeIntersection.Perp() > 0) { //suppress x,y=0,0 events
|
||||
if (PCQQQTimeCut) {
|
||||
plotter->Fill2D("PC_XY_Projection_QQQ_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, anodeIntersection.X(), anodeIntersection.Y(), "hPCQQQ");
|
||||
}
|
||||
plotter->Fill2D("PC_XY_Projection_QQQ" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100, anodeIntersection.X(), anodeIntersection.Y(), "hPCQQQ");
|
||||
}
|
||||
for (int j = i + 1; j < qqq.multi; j++)
|
||||
{
|
||||
if (qqq.id[i] == qqq.id[j])
|
||||
{
|
||||
int chWedge = -1;
|
||||
int chRing = -1;
|
||||
double eWedge = 0.0;
|
||||
double eWedgeMeV = 0.0;
|
||||
double eRing = 0.0;
|
||||
double eRingMeV = 0.0;
|
||||
double tRing = 0.0;
|
||||
int qqqID = -1;
|
||||
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
|
||||
{
|
||||
chWedge = qqq.ch[i];
|
||||
eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
|
||||
chRing = qqq.ch[j] - 16;
|
||||
eRing = qqq.e[j];
|
||||
tRing = static_cast<double>(qqq.t[j]);
|
||||
qqqID = qqq.id[i];
|
||||
}
|
||||
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
|
||||
{
|
||||
chWedge = qqq.ch[j];
|
||||
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
|
||||
chRing = qqq.ch[i] - 16;
|
||||
tRing = static_cast<double>(qqq.t[i]);
|
||||
eRing = qqq.e[i];
|
||||
qqqID = qqq.id[i];
|
||||
}
|
||||
else
|
||||
continue;
|
||||
|
||||
if (qqqCalibValid[qqq.id[i]][chRing][chWedge])
|
||||
{
|
||||
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
|
||||
eRingMeV = eRing * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
|
||||
}
|
||||
else
|
||||
continue;
|
||||
|
||||
// if (anodeIntersection.Z() != 0)
|
||||
{
|
||||
plotter->Fill2D("PC_Z_vs_QQQRing", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
|
||||
}
|
||||
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 2)
|
||||
{
|
||||
plotter->Fill2D("PC_Z_vs_QQQRing_2C", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
|
||||
plotter->Fill2D("PC_Z_vs_QQQRing_2C" + std::to_string(qqq.id[i]), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCzQQQ");
|
||||
plotter->Fill2D("PC_Z_vs_QQQWedge_2C", 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chWedge, "hPCzQQQ");
|
||||
}
|
||||
plotter->Fill2D("VertexRecon_QQQRingTC" + std::to_string(PCQQQTimeCut) + "PhiC" + std::to_string(PCQQQPhiCut), 600, -1300, 1300, 16, 0, 16, vector_closest_to_z.Z(), chRing, "hPCQQQ");
|
||||
double phi = TMath::ATan2(anodeIntersection.Y(), anodeIntersection.X()) * 180. / TMath::Pi();
|
||||
plotter->Fill2D("PolarAngle_Vs_QQQWedge" + std::to_string(qqqID), 360, -360, 360, 16, 0, 16, phi, chWedge, "hPCQQQ");
|
||||
// plotter->Fill2D("EdE_PC_vs_QQQ_timegate_ls1000"+std::to_string())
|
||||
|
||||
plotter->Fill2D("PC_Z_vs_QQQRing_Det" + std::to_string(qqqID), 600, -300, 300, 16, 0, 16, anodeIntersection.Z(), chRing, "hPCQQQ");
|
||||
//double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
|
||||
//double rho = 50. + 40. / 16. * (chRing + 0.5);
|
||||
|
||||
for (int k = 0; k < pc.multi; k++)
|
||||
{
|
||||
if(pc.index[k] >= 24)
|
||||
continue;
|
||||
|
||||
// double sinTheta = TMath::Sin((hitPos-vector_closest_to_z).Theta());
|
||||
double sinTheta = TMath::Sin((anodeIntersection-TVector3(0,0,90.0)).Theta());
|
||||
// double sinTheta = TMath::Sin((anodeIntersection-vector_closest_to_z).Theta());
|
||||
// double sinTheta = TMath::Sin((hitPos-TVector3(0,0,30.0)).Theta());
|
||||
// double sinTheta = TMath::Sin(hitPos.Theta());
|
||||
|
||||
if(cathodeHits.size()==2 && PCQQQPhiCut) {
|
||||
plotter->Fill2D("CalibratedQQQE_RvsCPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eRingMeV, pc.e[k]*sinTheta, "hPCQQQ");
|
||||
plotter->Fill2D("CalibratedQQQE_WvsCPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eWedgeMeV, pc.e[k]*sinTheta, "hPCQQQ");
|
||||
plotter->Fill2D("CalibratedQQQE_RvsPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eRingMeV, pc.e[k], "hPCQQQ");
|
||||
plotter->Fill2D("CalibratedQQQE_WvsPCE_TC" + std::to_string(PCQQQTimeCut) , 400, 0, 10, 400, 0, 30000, eWedgeMeV, pc.e[k], "hPCQQQ");
|
||||
plotter->Fill2D("PCQQQ_dTimevsdPhi", 200, -2000, 2000, 80, -200, 200, tRing - static_cast<double>(pc.t[k]), (hitPos.Phi()-anodeIntersection.Phi()) * 180. / TMath::Pi(), "hTiming");
|
||||
}
|
||||
|
||||
}
|
||||
}///qqq i==j case end
|
||||
} //j loop end
|
||||
} // qqq i loop end
|
||||
|
||||
TVector3 guessVertex(0,0,90.); //for run12, subtract anodeIntersection.Z() by ~74.0 seems to work
|
||||
//rho=40.0 mm is halfway between the cathodes(rho=42) and anodes(rho=37)
|
||||
double pcz_guess = 42.0/TMath::Tan((hitPos-guessVertex).Theta()) + guessVertex.Z(); //this is ideally kept to be all QQQ+userinput for calibration of pcz
|
||||
if(PCQQQTimeCut && PCQQQPhiCut && hitPos.Perp()>0 && anodeIntersection.Perp()>0 && cathodeHits.size()>=2) {
|
||||
plotter->Fill2D("pczguess_vs_qqqE",100,0,200,800,0,20,pcz_guess,qqqenergy,"pczguess");
|
||||
double pczoffset=30.0;
|
||||
//plotter->Fill2D("pczguess_vs_pcz_rad="+std::to_string(hitPos.Perp()),100,0,200,150,0,200,pcz_guess,anodeIntersection.Z(),"pczguess"); //entirely qqq-derived position vs entirely PC derived position
|
||||
plotter->Fill2D("pczguess_vs_pcz_phi="+std::to_string(hitPos.Phi()*180./M_PI),100,0,200,150,0,200,pcz_guess,anodeIntersection.Z()+pczoffset,"pczguess"); //entirely qqq-derived position vs entirely PC derived position
|
||||
plotter->Fill2D("pczguess_vs_pcz",100,0,200,150,0,200,pcz_guess,anodeIntersection.Z()+pczoffset);
|
||||
plotter->Fill2D("pcz_vs_pcPhi_rad="+std::to_string(hitPos.Perp()),360,0,360,150,0,200,anodeIntersection.Phi()*180./M_PI,anodeIntersection.Z()+pczoffset,"pczguess");
|
||||
}
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
// plotting sx3 strip hits vs anode phi
|
||||
if (sx3.ch[i] < 8 && anodeIntersection.Perp()>0)
|
||||
plotter->Fill2D("PCPhi_vs_SX3Strip", 100, -200, 200, 8 * 24, 0, 8 * 24, anodeIntersection.Phi() * 180. / TMath::Pi(), sx3.id[i] * 8 + sx3.ch[i]);
|
||||
}
|
||||
|
||||
if (anodeIntersection.Z() != 0 && cathodeHits.size() == 3)
|
||||
{
|
||||
plotter->Fill1D("PC_Z_proj_3C", 600, -300, 300, anodeIntersection.Z(), "hPCzQQQ");
|
||||
}
|
||||
|
||||
plotter->Fill2D("AnodeMaxE_Vs_Cathode_Sum_Energy", 2000, 0, 30000, 2000, 0, 30000, aEMax, cESum, "hGMPC");
|
||||
plotter->Fill1D("Correlated_Cathode_MaxAnode", 6, 0, 5, corrcatMax.size(), "hGMPC");
|
||||
plotter->Fill2D("Correlated_Cathode_VS_MaxAnodeEnergy", 6, 0, 5, 2000, 0, 30000, corrcatMax.size(), aEMax, "hGMPC");
|
||||
plotter->Fill1D("AnodeHits", 12, 0, 11, anodeHits.size(), "hGMPC");
|
||||
plotter->Fill2D("AnodeMaxE_vs_AnodeHits", 12, 0, 11, 2000, 0, 30000, anodeHits.size(), aEMax, "hGMPC");
|
||||
|
||||
if (anodeHits.size() < 1)
|
||||
{
|
||||
plotter->Fill1D("NoAnodeHits_CathodeHits", 6, 0, 5, cathodeHits.size(), "hGMPC");
|
||||
}
|
||||
|
||||
if(!det.valid && det.unmatched_front_chans.size()) {
|
||||
for(auto bch : det.valid_back_chans)
|
||||
for(auto fch : det.unmatched_front_chans) {
|
||||
if(det.frontR[fch].size() != 0) {
|
||||
double trial_left = det.back[bch].at(0) - det.frontR[fch].at(0);
|
||||
double frontR = det.frontR[fch].at(0);
|
||||
double backE = det.back[bch].at(0);
|
||||
double frontX = (trial_left - frontR)/(trial_left + frontR);
|
||||
plotter->Fill2D("bad_fronts_id"+std::to_string(id)+"_fR"+std::to_string(fch)+"_b"+std::to_string(bch),800,0,4096,800,0,4096,det.frontR[fch].at(0),det.back[bch].at(0),"bad_fronts_vb");
|
||||
plotter->Fill2D("be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(fch)+"R_b"+std::to_string(bch),200,-1,1,800,0,8192,
|
||||
frontX,backE,"evsx_bad_front");
|
||||
|
||||
//double backE = det.back[bch].at(0)*sx3BackGain[id][fch][bch];
|
||||
//double frontE = det.frontR[fch]
|
||||
}
|
||||
if(det.frontL[fch].size() != 0) {
|
||||
double trial_right = det.back[bch].at(0) - det.frontL[fch].at(0);
|
||||
double frontL = det.frontL[fch].at(0);
|
||||
double backE = det.back[bch].at(0);
|
||||
double frontX = (frontL - trial_right)/(frontL + trial_right);
|
||||
plotter->Fill2D("bad_fronts_id"+std::to_string(id)+"_fL"+std::to_string(fch)+"_b"+std::to_string(bch),800,0,4096,800,0,4096,det.frontL[fch].at(0),det.back[bch].at(0),"bad_fronts_vb");
|
||||
plotter->Fill2D("be_vs_x_sx3_id_"+std::to_string(id)+"_f"+std::to_string(fch)+"L_b"+std::to_string(bch),200,-1,1,800,0,8192,
|
||||
frontX,backE,"evsx_bad_front");
|
||||
}
|
||||
}//end for
|
||||
|
||||
}//end if(!det.valid)
|
||||
}//end for id
|
||||
}//end sx3.multi>1
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
void MakeVertexSX3::Terminate()
|
||||
{
|
||||
plotter->FlushToDisk();
|
||||
plotter->FlushToDisk(10);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -2,7 +2,7 @@
|
|||
export DATASET="27Al"
|
||||
export flip180="0"
|
||||
export flipa=0
|
||||
export anode_offset=1
|
||||
export anode_offset=0
|
||||
#declare -i run=28
|
||||
#while [[ $run -lt 34 ]]; do #runs 1 to 84
|
||||
# wrun=$(printf "%03d" $run)
|
||||
|
|
|
|||
18
run_sx3.sh
18
run_sx3.sh
|
|
@ -20,9 +20,9 @@ fi
|
|||
export DATASET="27Al"
|
||||
export flip180="0"
|
||||
#root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run09.root;
|
||||
if [[ 1 -eq 0 ]]; then
|
||||
if [[ 1 -eq 1 ]]; then
|
||||
#export timecut_low=230.0;
|
||||
#export timecut_low=400.0;
|
||||
export timecut_low=400.0;
|
||||
#export timecut_high=400.0;
|
||||
#unset timecut_low, timecut_high
|
||||
#export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run09.root;
|
||||
|
|
@ -30,17 +30,18 @@ if [[ 1 -eq 0 ]]; then
|
|||
#export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_011_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run11.root;
|
||||
export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_012_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run12.root;
|
||||
#export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_013_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run13.root;
|
||||
unset timecut_low
|
||||
#exit
|
||||
fi
|
||||
|
||||
#protons+gas, 27Al
|
||||
#export flip180="1"
|
||||
#export flip180="0"
|
||||
if [[ 1 -eq 1 ]] ; then
|
||||
export flipa=0
|
||||
export anode_offset=0
|
||||
export source_vertex=-200.0; #put the 'source' on the entrance window
|
||||
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_015_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run15.root;
|
||||
if [[ 1 -eq 0 ]] ; then
|
||||
#export flipa=0
|
||||
#export anode_offset=0
|
||||
#export source_vertex=-200.0; #put the 'source' on the entrance window
|
||||
#root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_015_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run15.root;
|
||||
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_017_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run17.root;
|
||||
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_018_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run18.root;
|
||||
root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_019_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run19.root;
|
||||
|
|
@ -58,7 +59,7 @@ fi
|
|||
#root -l -x results_run19.root results_run12.root -e "new TBrowser"
|
||||
#exit
|
||||
export DATASET="17F"
|
||||
export flip180="0"
|
||||
#export flip180="0"
|
||||
if [[ 1 -eq 0 ]]; then
|
||||
root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_005_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run05.root;
|
||||
root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_006_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run06.root;
|
||||
|
|
@ -82,6 +83,7 @@ export source_vertex=53.44; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/So
|
|||
export source_vertex=14.24; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_019_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run19.root;
|
||||
export source_vertex=-24.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_020_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run20.root;
|
||||
export source_vertex=-73.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run21.root;
|
||||
exit
|
||||
fi
|
||||
#17F reaction data
|
||||
#export flip180="0"
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
#include "testmodel.h"
|
||||
|
||||
int quit=0;
|
||||
void handler(int){quit=1;}
|
||||
void handler(int){quit=0;}
|
||||
|
||||
int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5};
|
||||
void scan_offset(){
|
||||
|
|
@ -26,18 +26,19 @@ void scan_offset(){
|
|||
auto c1=c.cd(1);
|
||||
c1->SetGrid(1,1);
|
||||
f = new TFile(Form("../../results_run%d.root",i));
|
||||
|
||||
if(i==12) {
|
||||
//TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int"));
|
||||
TH2F *h23 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int_A1C2"));
|
||||
if(h23) continue;
|
||||
std::cout << "aaa" << h23 << std::endl;
|
||||
h23->SetLineColorAlpha(kOrange,0.75);
|
||||
h23->Draw("box SAME");
|
||||
|
||||
h23->GetYaxis()->SetRangeUser(-200,200);
|
||||
h23->Draw("box");
|
||||
while(gPad->WaitPrimitive());
|
||||
} else {
|
||||
//TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int"));
|
||||
//TH2F *h2 = (TH2F*)(f->Get("pcz_vs_sx3pczguess_A1C2_strip12"));
|
||||
TH2F *h2 = (TH2F*)(f->Get("pcz_vs_sx3pczguess_A1C2"));
|
||||
std::cout << h2 << std::endl;
|
||||
//TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2"));
|
||||
if(!h2) continue;
|
||||
h2->SetTitle(Form("case%d",i));
|
||||
|
|
@ -48,6 +49,7 @@ void scan_offset(){
|
|||
}
|
||||
TF1 eqline("x","x",-200,200);
|
||||
eqline.Draw("SAME");
|
||||
|
||||
c1->Modified();
|
||||
c1->Update();
|
||||
ctr+=1;
|
||||
|
|
@ -68,13 +70,14 @@ void scan_offset(){
|
|||
while(gPad->WaitPrimitive());
|
||||
|
||||
files.emplace_back(f);
|
||||
std::cout <<"Test" << std::endl;
|
||||
if(i==21) {
|
||||
i=11;
|
||||
c.Clear();
|
||||
c.Divide(2,1);
|
||||
ctr=0;
|
||||
}
|
||||
if(quit) break;
|
||||
//if(quit) break;
|
||||
}
|
||||
for(auto file : files) {
|
||||
file->Close();
|
||||
|
|
|
|||
|
|
@ -31,7 +31,7 @@ void scan_offset_fix(){
|
|||
//TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int"));
|
||||
TH2F *h23 = (TH2F*)(f->Get("pczfix_vs_qqqpczguess_A1C2"));
|
||||
h23->SetLineColorAlpha(kOrange,0.75);
|
||||
h23->Draw("SAME");
|
||||
h23->Draw("box SAME");
|
||||
|
||||
} else {
|
||||
//TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int"));
|
||||
|
|
|
|||
|
|
@ -1,12 +1,12 @@
|
|||
#include <TF1.h>
|
||||
double model(double* x, double* p) {
|
||||
double result = x[0];
|
||||
double factor = 29.0;
|
||||
double factor = 0.0;
|
||||
double slope = 0.7;
|
||||
if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope;
|
||||
else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor;
|
||||
else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2;
|
||||
else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2;
|
||||
else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*3;
|
||||
return result;
|
||||
}
|
||||
|
||||
|
|
@ -17,7 +17,7 @@ double model_invert(double *y, double *q) {
|
|||
if(TMath::Abs(y[0]) < 16.2/slope) result = y[0]/slope;
|
||||
else if(TMath::Abs(y[0]) < 49.8/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor;
|
||||
else if(TMath::Abs(y[0]) < 85.2/slope ) result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2;
|
||||
else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*2;
|
||||
else result=y[0]/slope-TMath::Sign(1.0,y[0])*factor*3;
|
||||
return result+40;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -1 +0,0 @@
|
|||
pccal/pc_gm_coeffs.dat
|
||||
49
slope_intercept_results_17F.dat
Normal file
49
slope_intercept_results_17F.dat
Normal file
|
|
@ -0,0 +1,49 @@
|
|||
#Histogram Number Slope Intercept
|
||||
0 0.937314 -16.871
|
||||
1 0 0
|
||||
2 0.965461 -1.54376
|
||||
3 0.926501 -3.27662
|
||||
4 0.905634 2.54577
|
||||
5 0.905634 -11.0387
|
||||
6 0.853919 6.23079
|
||||
7 0.945588 -9.54044
|
||||
8 0.884454 -11.8262
|
||||
9 0.922501 -3.42538
|
||||
10 0.903053 9.28069
|
||||
11 0.914653 9.87642
|
||||
12 0.965332 13.2526
|
||||
13 0.923847 -3.41775
|
||||
14 0.93845 25.9901
|
||||
15 0.955424 12.324
|
||||
16 0.95116 4.99595
|
||||
17 0.910745 2.86648
|
||||
18 0.941376 4.57217
|
||||
19 0.871622 932.111
|
||||
20 1.00624 7.86358
|
||||
21 0.969834 -45.001
|
||||
22 0.89304 -31.5635
|
||||
23 0.933226 4.02193
|
||||
24 0 0
|
||||
25 0.941896 6.16135
|
||||
26 0.980284 2.86886
|
||||
27 0.983166 -3.82952
|
||||
28 0.978704 -2.89713
|
||||
29 0.964947 2.25786
|
||||
30 0.94514 0.925074
|
||||
31 0.977231 1.6493
|
||||
32 0.919527 5.82742
|
||||
33 0.972243 2.88061
|
||||
34 0.928892 7.61384
|
||||
35 0.947376 -0.644223
|
||||
36 0.875342 6.066
|
||||
37 0 0
|
||||
38 0.970953 6.262
|
||||
39 0 0
|
||||
40 0.918408 -3.27891
|
||||
41 0.913619 4.11288
|
||||
42 0.954083 2.21261
|
||||
43 0.993037 5.48924
|
||||
44 0 0
|
||||
45 0.926406 -19.719
|
||||
46 1.00459 5.14574
|
||||
47 0.942483 5.54183
|
||||
Loading…
Reference in New Issue
Block a user