From 2f43903269c554bbd15d056ca8b5056acfddbb11 Mon Sep 17 00:00:00 2001 From: Sudarsan B Date: Fri, 29 May 2026 17:45:57 -0400 Subject: [PATCH] * Main logic broken up into subunits that aim to be self-contained. The entirety of the analysis is aimed to follow: - find sx3 events, gainmatch them, calibrate, make a vector of these called SX3_Events - same step for QQQ, call it QQQ_Events - find anodeWire clusters, cathodeWire clusters, make vectors of these - make PC_Events from wire clusters, save 'anode only' and 'cathode only' cases just in case * since PC_Events, QQQ_Events, SX3_Events are all STL containers, we can pass them to functions that perform modular analyses * In this commit: - 4-wire offsets used, along with reliably figuring out phi-offset between QQQ/SX3/PC. We are close, except for some QQQ fine-tuning due to their angular extent not being 90 deg - the pczfix step now makes sense, 1-wire calculations also loosely match between guess and pcz - the nonlinearity correction/dynamic range fix is just adding to the resolution - p(a,a) data analysable by QQQ alphas show good kinematics, when doing the following gates: > A1C1 event pczs dithered, A1C2 events made into pczfix > phicut (45 deg) to gate out p+a correlations, SiE < 5 MeV && PCa > 6000 to select the alpha blob in p+a > Some selection on VertexReconXY so that the Perp() of the vertex is < 6mm > Selection on VertexReconZ so that z is in [-173.6, 100]. Fine alignment pending - Doing all the above gives reasonable p(a,a) kinematic curves with good statistics, Ex peaked at -0.7 MeV, close but not perfect - The above steps, when repeated with oneWire anode events stored in aClusters also yields a very reasonable kinematic locus, and Ex value, even more statistics - Not sure how much of this is autocorrelation stuff but > VertexReconZ vs Ef in QQQ (Ef is eloss-fixed alpha energy using path length) shows sensible trends. - Why Ex is not centered around zero might need more thought. - VertexReconXY is suspiciously well-centered, might need some more thought as well. * Some infrastructure that allows processing of 17F data is also in the pipeline now. * One fairly important bug got fixed which was ignoring qqq.id[0] when making QQQ_Events * Ideally, follow-ups to the above done on 27Al, 17F will make their own functions that are then called separately using booleans at the very top. * The fate of proton dE signals is out to jury. --- Armory/ClassDet.h | 4 +- Armory/ClassPW.h | 47 +- Armory/PC_StepLadder_Correction.h | 23 +- MakeVertex.C | 1333 ++++++++--------- MakeVertex.h | 5 - cutlist.txt | 4 + results/compare.C | 4 +- results/compareE.C | 1 - results/compareE_pcalpha.C | 26 + results/compareE_pcalpha_ow.C | 26 + results/haddnow | 6 + run_17F.sh | 56 +- run_27Al.sh | 27 +- run_sx3.sh | 32 +- scratch/sx3z_vs_phiz/compare_wire_offsets.C | 107 ++ .../sx3z_vs_phiz/compare_wire_offsets_c3.C | 88 ++ scratch/sx3z_vs_phiz/scan_offset.C | 89 +- scratch/sx3z_vs_phiz/scan_offset_a2c1.C | 86 ++ scratch/sx3z_vs_phiz/scan_offset_ainterp.C | 109 ++ scratch/sx3z_vs_phiz/scan_offset_fix.C | 9 +- scratch/sx3z_vs_phiz/scan_offset_fix_1d.C | 105 ++ .../sx3z_vs_phiz/scan_offset_fix_witha1c1.C | 88 ++ scratch/sx3z_vs_phiz/scan_offset_witha1c1.C | 87 ++ scratch/sx3z_vs_phiz/testmodel.h | 2 +- 24 files changed, 1517 insertions(+), 847 deletions(-) create mode 100644 cutlist.txt create mode 100644 results/compareE_pcalpha.C create mode 100644 results/compareE_pcalpha_ow.C create mode 100644 results/haddnow create mode 100644 scratch/sx3z_vs_phiz/compare_wire_offsets.C create mode 100644 scratch/sx3z_vs_phiz/compare_wire_offsets_c3.C create mode 100644 scratch/sx3z_vs_phiz/scan_offset_a2c1.C create mode 100644 scratch/sx3z_vs_phiz/scan_offset_ainterp.C create mode 100644 scratch/sx3z_vs_phiz/scan_offset_fix_1d.C create mode 100644 scratch/sx3z_vs_phiz/scan_offset_fix_witha1c1.C create mode 100644 scratch/sx3z_vs_phiz/scan_offset_witha1c1.C diff --git a/Armory/ClassDet.h b/Armory/ClassDet.h index 510cbb4..fa47e86 100644 --- a/Armory/ClassDet.h +++ b/Armory/ClassDet.h @@ -14,6 +14,7 @@ public: unsigned short ch[MAXMULTI]; unsigned short e[MAXMULTI]; unsigned long long t[MAXMULTI]; + unsigned long long tf[MAXMULTI]; unsigned short sn[MAXMULTI]; unsigned short digiCh[MAXMULTI]; @@ -28,6 +29,7 @@ public: ch[i] = 0; e[i] = 0; t[i] = 0; + tf[i] = 0; index[i] = 0; sn[i] = 0; digiCh[i] = 0; @@ -38,7 +40,7 @@ public: void Print(){ printf("=============================== multi : %u\n", multi); for( int i = 0; i < multi; i++) { - printf(" %3d | %2d-%-2d(%5d) %5u %15llu \n", i, id[i], ch[i], index[i], e[i], t[i]); + printf(" %3d | %2d-%-2d(%5d) %5u %15llu %15llu\n", i, id[i], ch[i], index[i], e[i], t[i], tf[i]); } } diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index 334b83e..41f12c4 100644 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -85,13 +85,14 @@ public: TVector3 GetAnodneMid(short id) const { return (An[id].first + An[id].second) * 0.5; } double GetAnodeTheta(short id) const { return (An[id].first - An[id].second).Theta(); } double GetAnodePhi(short id) const { return (An[id].first - An[id].second).Phi(); } + inline void PrintGeometry(); TVector3 GetCathodneMid(short id) const { return (Ca[id].first + Ca[id].second) * 0.5; } double GetCathodeTheta(short id) const { return (Ca[id].first - Ca[id].second).Theta(); } double GetCathodePhi(short id) const { return (Ca[id].first - Ca[id].second).Phi(); } void ClearHitInfo(); - void ConstructGeo(); + void ConstructGeo(int, int); void FindWireID(TVector3 pos, TVector3 direction, bool verbose = false); void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false); void CalTrack2(TVector3 sx3Pos, TVector3 anodeInt, bool verbose = false); @@ -144,7 +145,7 @@ inline void PW::ClearHitInfo() hitInfo.Clear(); } -inline void PW::ConstructGeo() +inline void PW::ConstructGeo(int aoffset=0, int coffset=0) { An.clear(); @@ -154,8 +155,8 @@ inline void PW::ConstructGeo() std::pair q1; // cathode double k = TMath::TwoPi()/24.; //48 solder thru holes, wires in every other one - double offset_a1 = -6*k-3*k; - double offset_c1 = -4*k -2*k - TMath::TwoPi()/48; //correct for a half-turn + double offset_a1 = -6*k-4*k;//-5*k-5*k; + double offset_c1 = -6*k+k/2.; //correct for a half-turn, 48 holes out of which 24 are used, anodes and cathodes staggered by half a turn. //-5*k-TMath::TwoPi()/48; //std::cerr << "Here!" << std::endl; //#include "../scratch/testing.h" double offset_a2 = offset_a1+wireShift*k; @@ -229,6 +230,36 @@ 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 void PW::PrintGeometry() { + for (size_t i = 0; i < An.size(); i++) + { + for (size_t j = 0; j < Ca.size(); j++) { + if ( Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190) { + continue; + } + + std::cout << i << ", " << j << ", " << Crossover[i][j][0].x << ", " << Crossover[i][j][0].y << ", " << Crossover[i][j][0].z << ", " << Crossover[i][j][1].x << ", " + << An[i].first.X() << ", " + << An[i].first.Y() << ", " + << An[i].first.Z() << ", " + + << Ca[j].first.X() << ", " + << Ca[j].first.Y() << ", " + << Ca[j].first.Z() << ", " + + << An[i].second.X() << ", " + << An[i].second.Y() << ", " + << An[i].second.Z() << ", " + + << Ca[j].second.X() << ", " + << Ca[j].second.Y() << ", " + << Ca[j].second.Z() + << std::endl; + } + std::cout << std::endl; + } +} + inline TVector3 PW::getClosestWirePosAtWirePhi(std::pair awire, double sx3phi_radian) { // 1. Get wire geometry TVector3 a1 = awire.first; // Top of the wire @@ -249,9 +280,13 @@ inline TVector3 PW::getClosestWirePosAtWirePhi(std::pair awi 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; + if(TMath::IsNaN(sx3phi_radian-test_pt.Phi())) { + //std::cout << "Moo" << std::endl; + 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 - + //std::cout << "Yee" << std::endl; + // If this is the smallest Delta Phi we've seen so far, save it! if (dPhi < min_delta_phi) { diff --git a/Armory/PC_StepLadder_Correction.h b/Armory/PC_StepLadder_Correction.h index 6618596..9ed25d9 100644 --- a/Armory/PC_StepLadder_Correction.h +++ b/Armory/PC_StepLadder_Correction.h @@ -8,23 +8,12 @@ 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(double *y, double *q) { - double result=y[0]; - double slope = 0.7; - 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*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}; + double slope = 0.56; + double z_grid[8] = {147.998,101.946,59.7634,19.6965,-19.6965,-59.7634,-101.946,-147.998}; //crossover z values 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 @@ -32,7 +21,7 @@ double model_invert(double* y, double* p) { break; } } - return result+80; + return result; } double model_a1c1(double* x, double* p) { @@ -48,12 +37,6 @@ double model_a1c1(double* x, double* p) { 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; } diff --git a/MakeVertex.C b/MakeVertex.C index 64163b7..e1c5687 100644 --- a/MakeVertex.C +++ b/MakeVertex.C @@ -27,6 +27,7 @@ Int_t colors[40] = { #include #include #include +#include #include #include @@ -43,6 +44,8 @@ double source_vertex = 53; //53 const double qqq_z = 100.0; const double anode_gain = 1.5146e-5; //channels --> MeV std::string dataset; +bool reactiondata=false; +bool Seven_MeV_Cut=false; TF1 pcfix_func("func",model_invert,-200,200); TF1 pcfix_func_a1c1("func_a1c1",model_invert_a1c1,-200,200); @@ -62,16 +65,14 @@ TGraph2D *qqqg=NULL, *crossoverg=NULL, *guessg=NULL; double z_to_crossover_rho(double z) { - return 9.20645e-5*z*z + 34.1973; -} +// return 9.20645e-5*z*z + 34.1973; + return 0.000165896*z*z + 4.61626e-08*z + 32.067; -double z_to_crossover_rho_cathode(double z) { - return 9.20645e-5*z*z + 34.1973; } // Global instances -PW pw_contr; PW pwinstance; +TRandom3 rnd_qqq, rnd_sx3; TVector3 hitPos; double qqqenergy, qqqtimestamp; class Event { @@ -118,7 +119,9 @@ bool qqqEcut; void protonAlphaHistograms(HistPlotter* plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events); - +void protonMiscHistograms(HistPlotter* plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events); +void protonMiscHistograms_sx3(HistPlotter* plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events); +void miscHistograms_oneWire(HistPlotter* plotter, std::vector QQQ_Events, std::vector>> aClusters); void MakeVertex::Begin(TTree * /*tree*/) { pcfix_func.SetNpx(100000); @@ -127,11 +130,7 @@ void MakeVertex::Begin(TTree * /*tree*/) plotter = new HistPlotter(option.Data(),"TFILE"); else plotter = new HistPlotter("Analyzer_SX3.root", "TFILE"); - - pw_contr.ConstructGeo(); - pwinstance.ConstructGeo(); - if(gROOT->IsBatch()) realtime=false; - + plotter->ReadCuts("cutlist.txt"); // --------------------------------------------------------- // 1. CRITICAL FIX: Initialize PC Arrays to Default (Raw) // --------------------------------------------------------- @@ -140,7 +139,8 @@ void MakeVertex::Begin(TTree * /*tree*/) pcSlope[i] = 1.0; // Default slope = 1 (preserves Raw energy) pcIntercept[i] = 0.0; // Default intercept = 0 } - + rnd_qqq.SetSeed(0); + rnd_sx3.SetSeed(0); if(getenv("DATASET")) dataset = std::string(getenv("DATASET")); if(getenv("source_vertex")) @@ -148,15 +148,22 @@ void MakeVertex::Begin(TTree * /*tree*/) std::cout << "Dataset set to " << dataset << std::endl; std::cout << "source_vertex set to " << source_vertex << std::endl; + int aoffset = 0; + int coffset = 0; + if(getenv("anode_offset")) { + aoffset = std::atoi(getenv("anode_offset")); + std::cout << "Offseting anodes by " << aoffset << " wires." << std::endl; + } + if(getenv("cathode_offset")) { + coffset = std::atoi(getenv("cathode_offset")); + std::cout << "Offseting cathodes by " << coffset << " wires." << std::endl; + } + pwinstance.ConstructGeo(aoffset,coffset); + //pwinstance.PrintGeometry(); - if(getenv("flipa")) { - int flip_offset = std::atoi(getenv("anode_offset")); - int yes_to_flip = std::atoi(getenv("flipa")); - if(yes_to_flip && flip_offset) { - std::cout << "Flipping anodes and offseting by " << flip_offset << " wires." << std::endl; - } else if(flip_offset){ - std::cout << "Offseting anodes without flip by " << flip_offset << " wires." << std::endl; - } + if(getenv("reactiondata")) { + reactiondata = std::atoi(getenv("reactiondata")); + std::cout << "Analyzing dataset as reactiondata" << std::endl; } fflush(stdout); @@ -251,62 +258,14 @@ void MakeVertex::Begin(TTree * /*tree*/) //cm_to_MeV.Eval(MeV_to_cm.Eval(detectedE)-PathLength) gives energy of particle before it traversed 'path length' - if(realtime) { - can1 = new TCanvas("wireindex","c1",0,0,640,480); - can2 = new TCanvas("3d","c2",650,0,640,480); - can1->cd(); - //can2->SetFillColor(30); - frame = new TH3D("frame","frame",1000,-100,100,1000,-100,100,1000,-200,200); - hha =new TH1F("hha","Anode Ecal vs wire#",48,-12,36); - hhc =new TH1F("hhc","Cathode Ecal vs wire#",48,-12,36); - hha->SetLineColor(kRed); - hha->GetYaxis()->SetRangeUser(0,16384); - hha->GetXaxis()->SetTitle("press any key, interrupt/refresh or double click to continue.."); - hha->Draw(); - hhc->Draw("SAME"); - can1->Modified(); - can1->Update(); - can1->BuildLegend(); - can2->cd(); - frame->Draw(); - for(int i=0; i<24; i++) { - plc[i] = new TPolyLine3D(2); - pla[i] = new TPolyLine3D(2); - pla[i]->SetPoint(0,pwinstance.An[i].first.X(),pwinstance.An[i].first.Y(),pwinstance.An[i].first.Z()); - pla[i]->SetPoint(1,pwinstance.An[i].second.X(),pwinstance.An[i].second.Y(),pwinstance.An[i].second.Z()); - plc[i]->SetPoint(0,pwinstance.Ca[i].first.X(),pwinstance.Ca[i].first.Y(),pwinstance.Ca[i].first.Z()); - plc[i]->SetPoint(1,pwinstance.Ca[i].second.X(),pwinstance.Ca[i].second.Y(),pwinstance.Ca[i].second.Z()); - plc[i]->SetLineStyle(kDotted); - pla[i]->SetLineStyle(kDotted); - pla[i]->SetLineWidth(1.); - plc[i]->SetLineWidth(1.); - plc[i]->Draw("same"); - pla[i]->Draw("same"); - plc[i]->SetLineColor(colors[i]); - pla[i]->SetLineColor(colors[i]); - } - crossoverg = new TGraph2D(1); - crossoverg->SetName("crossoverg"); - crossoverg->SetMarkerStyle(20); - crossoverg->SetMarkerColor(kBlue+3); - qqqg = new TGraph2D(1); - qqqg->SetName("qqqg"); - qqqg->SetMarkerColor(kRed); - qqqg->SetMarkerStyle(42); - - crossoverg->SetPoint(0,0,0,0); - qqqg->SetPoint(0,0,0,qqq_z); - crossoverg->Draw("P same"); - qqqg->Draw("P same"); - - trajectory=new TPolyLine3D(2); - trajectory->SetPoint(0,0,0,0); - trajectory->SetPoint(1,0,0,0); - trajectory->Draw("same"); - - can2->Modified(); - can2->Update(); - } + if(dataset=="17F" && reactiondata) { + fChain->SetBranchAddress("miscMulti", &misc.multi, &b_miscMulti); + fChain->SetBranchAddress("miscID", &misc.id, &b_miscID); + fChain->SetBranchAddress("miscCh", &misc.ch, &b_miscCh); + fChain->SetBranchAddress("miscE", &misc.e, &b_miscE); + fChain->SetBranchAddress("miscT", &misc.t, &b_miscT); + fChain->SetBranchAddress("miscF", &misc.tf, &b_miscTf); + } } Bool_t MakeVertex::Process(Long64_t entry) @@ -331,7 +290,15 @@ Bool_t MakeVertex::Process(Long64_t entry) b_pcCh->GetEntry(entry); b_pcE->GetEntry(entry); b_pcT->GetEntry(entry); - + if(dataset=="17F" && reactiondata) { + b_miscMulti->GetEntry(entry); + b_miscID->GetEntry(entry); + b_miscCh->GetEntry(entry); + b_miscE->GetEntry(entry); + b_miscT->GetEntry(entry); + b_miscTf->GetEntry(entry); + } + double timecut_low = getenv("timecut_low")?std::atof(getenv("timecut_low")):0; double timecut_high = getenv("timecut_high")?std::atof(getenv("timecut_high")):1e15 ; @@ -368,7 +335,7 @@ Bool_t MakeVertex::Process(Long64_t entry) double value=sx3.e[i]; int gch = sx3.id[i]*4+(sx3.ch[i]-8); if(id<12) Fsx3.at(id).fillevent("BACK",sx3ch,value); - Fsx3.at(id).ts = static_cast(sx3.t[i]); + Fsx3.at(id).ts = static_cast(sx3.t[i])+(rnd_sx3.Uniform(16.0)-8.0); plotter->Fill2D("sx3backs_all_raw",100,0,100,800,0,4096,gch,sx3.e[i]); } else { int sx3ch=sx3.ch[i]/2; @@ -411,29 +378,19 @@ Bool_t MakeVertex::Process(Long64_t entry) // Note that this will be different for the upstream barrel, when it gets implemented double backE = det.backE*sx3BackGain[id][det.stripF][det.stripB]; //if(backE<2000) continue; - det.stripF=3-det.stripF; - - 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 + //det.stripF=3-det.stripF; + if(Seven_MeV_Cut && backE<5000) continue; + double alpha_n = TMath::ATan2((2*(3-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") { - //if(dataset=="17F") - // phi_n+=180;//run 37 in 17F--> - //} - //} + 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 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(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); - - //plotter->Fill2D("SX3CartesianPlot", 200, -100, 100, 200, -100, 100, 88.0*TMath::Cos(phi_n),88.0*TMath::Sin(phi_n), "hCalSX3"); plotter->Fill2D("SX3CartesianPlot" + std::to_string(id), 200, -100, 100, 200, -100, 100, 88.0*TMath::Cos(phi_n),88.0*TMath::Sin(phi_n), "hCalSX3"); - } if(det.valid && det.stripF!=DEFAULT_NULL && det.stripB!=DEFAULT_NULL) { plotter->Fill2D("sx3backs_raw",100,0,100,800,0,8192,det.stripB+4*id,det.backE); @@ -509,16 +466,16 @@ Bool_t MakeVertex::Process(Long64_t entry) 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(qqq.t[j]); - tWedge = static_cast(qqq.t[i]); + tRing = static_cast(qqq.t[j])+(rnd_qqq.Uniform(16.0)-8.0); + tWedge = static_cast(qqq.t[i])+(rnd_qqq.Uniform(16.0)-8.0); } 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(qqq.t[i]); - tWedge = static_cast(qqq.t[j]); + tRing = static_cast(qqq.t[i])+(rnd_qqq.Uniform(16.0)-8.0); + tWedge = static_cast(qqq.t[j])+(rnd_qqq.Uniform(16.0)-8.0); } else continue; @@ -534,18 +491,21 @@ Bool_t MakeVertex::Process(Long64_t entry) if(eRingMeV/eWedgeMeV > 3.0 || eRingMeV/eWedgeMeV<1.0/3.0) continue; //if(eRingMeV<1.2 || eWedgeMeV<1.2) continue; - - double theta = 2 * TMath::Pi() * (-qqq.id[i] * 16 + (15-chWedge) + 0.5)/(16*4); + + //double theta = 2 * TMath::Pi() * (-qqq.id[i] * 16 + (15-chWedge) + 0.5)/(16*4); + double theta = (M_PI/180.)*(-90*qqq.id[i]+(87./16.)*((15-chWedge)+0.5)+3.0); 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),qqq_z), 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),qqq_z), eRing, eWedge, tRing, tWedge,chRing+qqq.id[i]*16, chWedge+qqq.id[i]*16); - if(qqq.id[i]>=1) { - QQQ_Events.push_back(qqqevent); - QQQ_Events_Raw.push_back(qqqeventr); - plotter->Fill2D("WedgeE_Vs_RingECal_selected", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ"); - } + assert(qqq.id[i]>=0); + if(Seven_MeV_Cut &&(eRingMeV<6.6 || eWedgeMeV< 6.6)) continue; + QQQ_Events.push_back(qqqevent); + QQQ_Events_Raw.push_back(qqqeventr); + plotter->Fill2D("WedgeE_Vs_RingECal_selected", 1000, 0, 10, 1000, 0, 10, eWedgeMeV, eRingMeV, "hCalQQQ"); 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"); 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"); @@ -585,13 +545,13 @@ Bool_t MakeVertex::Process(Long64_t entry) 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 theta = 2 * TMath::Pi() * (-qqq.id[i] * 16 + (15-chWedge) + 0.5)/(16*4); + //double theta = 2 * TMath::Pi() * (-qqq.id[i] * 16 + (15-chWedge) + 0.5)/(16*4); + double theta = -90.+90*qqq.id[i]+(87./16.)*((15-chWedge)+0.5)+3.0; + theta *= (M_PI/180.); double rho = 50. + (50. / 16.) * (chRing + 0.5); //"?" double x = rho * TMath::Cos(theta); double y = rho * TMath::Sin(theta); hitPos.SetXYZ(x, y, qqq_z); - //if(realtime) qqqg->SetPoint(0,hitPos.X(),hitPos.Y(),hitPos.Z()); - if(realtime) qqqg->AddPoint(hitPos.X(),hitPos.Y(),hitPos.Z()); qqqenergy = eRingMeV; qqqtimestamp = tRing; HitNonZero = true; @@ -605,13 +565,12 @@ Bool_t MakeVertex::Process(Long64_t entry) typedef std::unordered_map> WireEvent; //this stores nearest neighbour wire events, or a 'cluster' WireEvent aWireEvents, cWireEvents; //naming for book keeping - + aWireEvents.clear(); aWireEvents.reserve(24); - if(realtime) { - hha->Reset(); - hhc->Reset(); - } + cWireEvents.clear(); + cWireEvents.reserve(24); + // PC Gain Matching and Filling double anodeT = -99999; double cathodeT = 99999; @@ -620,7 +579,7 @@ Bool_t MakeVertex::Process(Long64_t entry) for (int i = 0; i < pc.multi; i++) { //std::cout << pc.index[i] << " " << pc.e[i] << " " << std::endl; - if (pc.e[i] > 10) + if (pc.e[i] > 25) { plotter->Fill2D("PC_Index_Vs_Energy", 48, 0, 48, 2000, 0, 30000, pc.index[i], static_cast(pc.e[i]), "hRawPC"); } else { @@ -637,37 +596,13 @@ Bool_t MakeVertex::Process(Long64_t entry) { anodeT = static_cast(pc.t[i]); anodeIndex = pc.index[i]; - - if(getenv("flipa")) { - int flip_offset = std::atoi(getenv("anode_offset")); - int yes_to_flip = std::atoi(getenv("flipa")); - if(yes_to_flip && flip_offset) { - int flipped_index = (23-anodeIndex+flip_offset)%24; - aWireEvents[flipped_index] = std::tuple(flipped_index,pc.e[i],static_cast(pc.t[i])); - //std::cout << "Flipping anodes and offseting by " << flip_offset << " wires." << std::endl; - } else if(flip_offset){ - int offset_index = (anodeIndex+flip_offset)%24; - aWireEvents[pc.index[i]] = std::tuple(offset_index,pc.e[i],static_cast(pc.t[i])); - //std::cout << "Offseting anodes without flip by " << offset_index << " wires." << std::endl; - } else - aWireEvents[pc.index[i]] = std::tuple(pc.index[i],pc.e[i],static_cast(pc.t[i])); - } else - aWireEvents[pc.index[i]] = std::tuple(pc.index[i],pc.e[i],static_cast(pc.t[i])); - if(realtime) hha->SetBinContent(hha->FindFixBin(anodeIndex),pc.e[i]); + aWireEvents[pc.index[i]] = std::tuple(pc.index[i],pc.e[i],static_cast(pc.t[i])); } else { cathodeT = static_cast(pc.t[i]); cathodeIndex = pc.index[i] - 24; - if(getenv("flipc")) { - int flip_offset = std::atoi(getenv("flipc")); - int flipped_index = (cathodeIndex+flip_offset)%24; - cWireEvents[flipped_index] = std::tuple(flipped_index,pc.e[i],static_cast(pc.t[i])); - } else { - cWireEvents[pc.index[i]-24] = std::tuple(pc.index[i]-24,pc.e[i],static_cast(pc.t[i])); - } - - if(realtime) hhc->SetBinContent(hhc->FindFixBin(cathodeIndex),pc.e[i]); + cWireEvents[pc.index[i]-24] = std::tuple(pc.index[i]-24,pc.e[i],static_cast(pc.t[i])); } if (anodeT != -99999 && cathodeT != 99999) @@ -742,11 +677,13 @@ Bool_t MakeVertex::Process(Long64_t entry) std::vector>> cClusters = pwinstance.Make_Clusters(cWireEvents); std::vector> sumE_AC; + int actr = 0; for(auto aCluster: aClusters) { + int cctr = 0; for(auto cCluster: cClusters) { //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 && aCluster.size() && cCluster.size()) { + if(alpha!=9999999 && apSumE!=-1 && aCluster.size() && cCluster.size()) { //needs both anodes and cathodes, AND for the crossover to fall in [-173.6, 173.6] //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. @@ -755,34 +692,81 @@ 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 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. + } + if(cCluster.size()!=0 && actr==0) { + Event PCEvent_OnlyC(crossover,-1,cpMaxE,-1,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; - } + } + + if(aCluster.size()!=0 && cctr==0) { //avoid double-counting + Event PCEvent_OnlyA(crossover,apSumE,-1,apTSMaxE,-1); //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); + } + cctr++; } + actr++; } - - + + bool is_oxygen=false; + TRandom3 rnd; + rnd.SetSeed();//random seed set + if(dataset=="17F" && reactiondata) { + int ctr=0; + for(auto qqqevent: QQQ_Events) { + double ts_rf = -987654321; + double ts_needle = -987654321; + double ts_mcp =-987654321; + double ts_qqq = static_cast(qqqevent.Time1) + (rnd.Uniform(16.0)-8.0); + bool found_rf=false; + bool found_mcp=false; + bool found_needle=false; + for(int j=0; j< misc.multi; j++) { + plotter->Fill1D("channels_misc",20,0,20,misc.ch[j],"misc"); + if(misc.ch[j]==2) { //Needle + plotter->Fill2D("needle_vs_qqqE",800,0,16384,800,0,10,misc.e[j],qqqevent.Energy1,"misc"); + ts_needle = static_cast(misc.t[j])+static_cast(misc.tf[j]); + found_needle=1; + plotter->Fill1D("dt_qqq_needle",800,-2000,2000,ts_qqq-ts_needle,"misc"); + } + if(misc.ch[j]==3) { //RF + ts_rf = static_cast(misc.t[j])+static_cast(misc.tf[j]); + found_rf=1; + plotter->Fill1D("dt_qqq_rf",800,-2000,2000,ts_qqq-ts_rf,"misc"); + } + if(misc.ch[j]==4) { //mcp + ts_mcp = static_cast(misc.t[j])+static_cast(misc.tf[j]); + found_mcp=1; + plotter->Fill1D("dt_qqq_mcp",800,-2000,2000,ts_qqq-ts_mcp,"misc"); + } + } + if(found_rf && found_mcp) { + if(ctr==0) plotter->Fill1D("dt_rf_mcp",500,-1000,1000,ts_rf-ts_mcp,"misc"); + double dt_rf_mcp = ts_rf - ts_mcp; + double dt_qqq_rf = ts_qqq-ts_rf; + double dt_qqq_mcp = ts_qqq-ts_mcp; + plotter->Fill2D("dt(qqq,rf)_vs_(rf,mcp)",640,-2000,2000,640,-2000,2000,dt_qqq_rf,dt_rf_mcp,"misc"); + plotter->Fill2D("dt_(qqq,mcp)_vs_(qqq,rf)",640,-1400,2000,640,-2000,2000,dt_qqq_mcp,dt_qqq_rf,"misc"); + plotter->Fill2D("dt_(qqq,mcp)_vs_(rf,mcp)",640,-1400,-600,640,-2000,2000,dt_qqq_mcp,dt_rf_mcp,"misc"); + + } + ctr+=1; + } + } + if(process_alpha_proton_scattering) { protonAlphaHistograms(plotter,QQQ_Events,SX3_Events,PC_Events); //return kTRUE; }//end if(process_alpha_proton_scattering) + protonMiscHistograms(plotter,QQQ_Events,SX3_Events,PC_Events); + //protonMiscHistograms_sx3(plotter,QQQ_Events,SX3_Events,PC_Events); + miscHistograms_oneWire(plotter, QQQ_Events, aClusters); + return kTRUE; - - if(QQQ_Events.size() && PC_Events.size()) + if(QQQ_Events.size() && PC_Events.size()) plotter->Fill2D("PCEv_vs_QQQEv",20,0,20,20,0,20,QQQ_Events.size(),PC_Events.size()); plotter->Fill2D("ac_vs_cc",20,0,20,20,0,20,aClusters.size(),cClusters.size(),"wiremult"); @@ -799,89 +783,111 @@ Bool_t MakeVertex::Process(Long64_t entry) for(auto sx3event: SX3_Events) { + /*for(const auto pcevent : PC_Events_OnlyAnode) { + plotter->Fill1D("dt_oa_sx3",800,-2000,2000,sx3event.Time1 - apTSMaxE,"pcev_onlyanode"); + if(sx3event.Time1 - apTSMaxE < 150) { + plotter->Fill1D("dt_oa_sx3_gated",800,-2000,2000,sx3event.Time1 - pcevent.Energy1,"pcev_onlyanode"); + plotter->Fill2D("dt_oa_sx3_gated_vs_sx3E",800,-2000,2000,800,0,10,sx3event.Time1 - pcevent.Time1,sx3event.Energy1,"pcev_onlyanode"); + plotter->Fill2D("dEa_oa_Esx3_TC1_ignC"+std::to_string(acluster.size()),400,0,10,800,0,40000,sx3event.Energy1,pcevent.Energy1,"pcev_onlyanode"); + plotter->Fill2D("pcPhi_oa_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,"pcev_onlyanode"); + plotter->Fill2D("pcZ_oa_sx3Z_TC1_ignC"+std::to_string(acluster.size()),300,-100,200,400,-200,200,sx3event.pos.Z(),pc_closest.Z(),"pcev_onlyanode"); + + double sx3theta = TMath::ATan2(88.0,sx3event.pos.Z()-source_vertex); + double pczguess = z_to_crossover_rho(pc_closest.Z())/TMath::Tan(sx3theta) + source_vertex; + double sinTheta = TMath::Sin(sx3theta); + + plotter->Fill2D("pcZ_oa_sx3pczguess_TC1_ignC"+std::to_string(acluster.size()),300,-100,200,400,-200,200,pczguess,pc_closest.Z(),"pcev_onlyanode"); + 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; + + double sinTheta2 = TMath::Sin(TMath::ATan2(88.0,sx3event.pos.Z()-vector_closest_to_axis.Z())); + plotter->Fill2D("dEa3_oa_Esx3_TC1_ignC"+std::to_string(acluster.size()),400,0,10,800,0,30000,sx3event.Energy1,apSumE*sinTheta2,"pcev_ + + + "); + + plotter->Fill2D("vertexZ_oa_sx3Z_TC1_ignC"+std::to_string(acluster.size()),300,-100,200,400,-200,200,sx3event.pos.Z(),vector_closest_to_axis.Z(),"pcev_onlyanode"); + plotter->Fill2D("vertexXY_oa_TC1_ignC"+std::to_string(acluster.size()),200,-100,100,200,-100,100,vector_closest_to_axis.X(),vector_closest_to_axis.Y(),"pcev_onlyanode"); + } + }*/ + if(sx3event.Energy1<1.2) continue; 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"); - + if(sx3event.Time1 - apTSMaxE < 150) { + bool phicut = sx3event.pos.Phi() <= pc_closest.Phi()+TMath::Pi()/4. && sx3event.pos.Phi() >= pc_closest.Phi()-TMath::Pi()/4.; double sx3theta = TMath::ATan2(88.0,sx3event.pos.Z()-source_vertex); - double pczguess = 37.0/TMath::Tan(sx3theta) + source_vertex; + double pczguess = z_to_crossover_rho(pc_closest.Z())/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"); - + + plotter->Fill1D("dt_anode_ainterp_sx3_gated",800,-2000,2000,sx3event.Time1 - apTSMaxE,"ainterp_noc"); + plotter->Fill2D("dt_anode_ainterp_sx3_gated_vs_sx3E",800,-2000,2000,800,0,10,sx3event.Time1 - apTSMaxE,sx3event.Energy1,"ainterp_noc"); + plotter->Fill2D("dEa_ainterp_Esx3_TC1_ignC_a"+std::to_string(acluster.size()),400,0,10,800,0,40000,sx3event.Energy1,apSumE,"ainterp_noc"); + plotter->Fill2D("pcPhi_ainterp_sx3Phi_TC1_ignC_a"+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_a"+std::to_string(acluster.size())+"_PC"+std::to_string(phicut),300,-100,200,400,-200,200,sx3event.pos.Z(),pc_closest.Z(),"ainterp_noc"); + + + plotter->Fill2D("pcZ_ainterp_sx3pczguess_TC1_ignC_a"+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"); + double sinTheta2 = TMath::Sin((sx3event.pos - vector_closest_to_axis).Theta());; + plotter->Fill2D("dEa3_ainterp_Esx3_TC1_ignC_a"+std::to_string(acluster.size())+"_PC"+std::to_string(phicut),1200,0,30,800,0,30000,sx3event.Energy1,apSumE*sinTheta2*3.,"ainterp_noc"); + + plotter->Fill2D("vertexZ_ainterp_sx3Z_TC1_ignC_a"+std::to_string(acluster.size()),300,-100,200,800,-400,400,sx3event.pos.Z(),vector_closest_to_axis.Z(),"ainterp_noc"); + plotter->Fill1D("vertexZ1d_ainterp_sx3Z_TC1_ignC_a"+std::to_string(acluster.size()),800,-400,400,vector_closest_to_axis.Z(),"ainterp_noc"); + plotter->Fill2D("vertexXY_ainterp_TC1_ignC_a"+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]; - if(sx3event.Time1 -(double)std::get<2>(awire)< -150) { + if(sx3event.Time1 -(double)std::get<2>(awire)< 150) { //plotter->Fill2D("sx3_z_phi2_awire"+std::to_string(std::get<0>(awire)), 400,-100,100, 100, -200,200,sx3event.pos.Z(), sx3event.pos.Phi()*180/M_PI ); //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("anodeNum_vs_stripNum",64,0,64,24,0,24,sx3event.ch2,i,"onewire"); + bool sx3diagonalphi = (!plotter->FindCut("anode_sx3_diag1")->IsInside(sx3event.ch2,i)) || (!plotter->FindCut("anode_sx3_diag2")->IsInside(sx3event.ch2,i)); + plotter->Fill2D("anodeNum_vs_stripNum_diag"+std::to_string(sx3diagonalphi),64,0,64,24,0,24,sx3event.ch2,i,"onewire"); + plotter->Fill2D("onewire_dEa_Esx3_TC1_fullev"+std::to_string(PC_Events.size()>0)+"_PC"+std::to_string(sx3diagonalphi),400,0,10,800,0,40000,sx3event.Energy1,std::get<1>(awire),"onewire"); + //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,"onewire"); //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 ); } } - + if(cWireEvents.find(i) != cWireEvents.end()) { auto cwire = cWireEvents[i]; - if(sx3event.Time1 -(double)std::get<2>(cwire) < -150) { + if(sx3event.Time1 -(double)std::get<2>(cwire) < 150) { //plotter->Fill2D("sx3_z_phi2_cwire"+std::to_string(std::get<0>(cwire)),400,-100,100, 100, -200,200,sx3event.pos.Z(), sx3event.pos.Phi()*180/M_PI ); //plotter->Fill2D("sx3_z_strip#_cwire"+std::to_string(std::get<0>(cwire)),400,-100,100, 100, -50,50,sx3event.pos.Z(), sx3event.ch2 ); - plotter->Fill2D("onewire_dEc_Esx3_fullev"+std::to_string(PC_Events.size()>0),400,0,10,800,0,40000,sx3event.Energy1,std::get<1>(cwire)); - plotter->Fill2D("onewire_cNum_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("onewire_dEc_Esx3_fullev"+std::to_string(PC_Events.size()>0),400,0,10,800,0,40000,sx3event.Energy1,std::get<1>(cwire),"onewire"); + plotter->Fill2D("onewire_cNum_sx3Phi_TC1_fullev"+std::to_string(PC_Events.size()>0),24,0,24,120,-360,360,i,sx3event.pos.Phi()*180./M_PI,"onewire"); } } }//for 'i' loop } - - for(auto qqqevent: QQQ_Events) { - for(int i=0; i<24; i++) { - if(aWireEvents.find(i) != aWireEvents.end()) { - auto awire = aWireEvents[i]; - if(qqqevent.Time1 -(double)std::get<2>(awire)< -150) { - plotter->Fill2D("onewire_dEa_Eqqq_TC1_fullev"+std::to_string(PC_Events.size()>0),400,0,10,800,0,40000,qqqevent.Energy1,std::get<1>(awire)); - plotter->Fill2D("onewire_aNum_QQQPhi_TC1_fullev"+std::to_string(PC_Events.size()>0),24,0,24,120,-360,360,i,qqqevent.pos.Phi()*180./M_PI); - } - } - - if(cWireEvents.find(i) != cWireEvents.end()) { - auto cwire = cWireEvents[i]; - if(qqqevent.Time1 -(double)std::get<2>(cwire) < -150) { - plotter->Fill2D("onewire_dEc_Eqqq_TC1_fullev"+std::to_string(PC_Events.size()>0),400,0,10,800,0,40000,qqqevent.Energy1,std::get<1>(cwire)); - plotter->Fill2D("onewire_cNum_QQQPhi_TC1_fullev"+std::to_string(PC_Events.size()>0),24,0,24,120,-360,360,i,qqqevent.pos.Phi()*180./M_PI); - } - } - }//for 'i' loop - } - + + + + + for(auto sx3event:SX3_Events) { PCSX3TimeCut=false; 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 + if(sx3event.Time1 - pcevent.Time1 < 200)//-150 for alphas PCASX3TimeCut = 1; - if(sx3event.Time1 - pcevent.Time2 < 0)//-200 for alphas + if(sx3event.Time1 - pcevent.Time2 < 200)//-200 for alphas PCCSX3TimeCut = 1; PCSX3TimeCut = PCASX3TimeCut && PCCSX3TimeCut; - + bool phicut = sx3event.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/4. && sx3event.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/4.; plotter->Fill1D("dt_pcA_sx3B",640,-2000,2000,sx3event.Time1 - pcevent.Time1); @@ -893,7 +899,8 @@ Bool_t MakeVertex::Process(Long64_t entry) if(pcevent.multi1==1 && pcevent.multi2==2) plotter->Fill2D("dE_E_Cathodesx3B_a1c2",400,0,30,800,0,10000,sx3event.Energy1,pcevent.Energy2); if(pcevent.multi1==2 && pcevent.multi2==1) plotter->Fill2D("dE_E_Anodesx3B_a2c1",400,0,30,800,0,40000,sx3event.Energy1,pcevent.Energy1); if(pcevent.multi1==2 && pcevent.multi2==1) plotter->Fill2D("dE_E_Cathodesx3B_a2c1",400,0,30,800,0,10000,sx3event.Energy1,pcevent.Energy2); - plotter->Fill2D("sx3phi_vs_pcphi"+std::to_string(sx3event.Time1 - pcevent.Time1<-150),100,-360,360,100,-360,360,sx3event.pos.Phi()*180/M_PI,pcevent.pos.Phi()*180/M_PI); + plotter->Fill2D("sx3phi_vs_pcphi"+std::to_string(sx3event.Time1 - pcevent.Time1<150),100,-360,360,100,-360,360,sx3event.pos.Phi()*180/M_PI,pcevent.pos.Phi()*180/M_PI); + plotter->Fill1D("d_sx3phi_minus_pcphi"+std::to_string(sx3event.Time1 - pcevent.Time1<150),180,-360,360,sx3event.pos.Phi()*180/M_PI-pcevent.pos.Phi()*180/M_PI); if(PCSX3TimeCut) { plotter->Fill1D("dt_pcA_sx3B_timecut",640,-2000,2000,sx3event.Time1 - pcevent.Time1); plotter->Fill1D("dt_pcC_sx3B_timecut",640,-2000,2000,sx3event.Time1 - pcevent.Time2); @@ -901,26 +908,37 @@ Bool_t MakeVertex::Process(Long64_t entry) plotter->Fill2D("xyplot_sx3"+std::to_string(sx3event.ch2/4),100,-100,100,100,-100,100,pcevent.pos.X(),pcevent.pos.Y()); plotter->Fill2D("pcz_vs_pcphi_TimeCut",600,-200,200,120,-360,360,pcevent.pos.Z(),pcevent.pos.Phi()*180/M_PI); //x-axis is all Si det, y-axis is PC anode+cathode only } - + double sx3rho = 88.0;//approximate barrel radius double sx3z = sx3event.pos.Z(); //w.r.t target origin at 90 for run12 double pcz = pcevent.pos.Z(); double calcsx3theta = TMath::ATan2(sx3rho-z_to_crossover_rho(pcz),sx3z-pcz); - plotter->Fill2D("dE2_E_Anodesx3B",400,0,30,800,0,40000,sx3event.Energy1,pcevent.Energy1*TMath::Sin(calcsx3theta)); - plotter->Fill2D("dE2_E_Cathodesx3B",400,0,30,800,0,10000,sx3event.Energy1,pcevent.Energy2*TMath::Sin(calcsx3theta)); - - + plotter->Fill2D("dE4_E_Anodesx3B",400,0,30,800,0,40000,sx3event.Energy1,pcevent.Energy1*TMath::Sin(calcsx3theta)); + plotter->Fill2D("dE4_E_Cathodesx3B",400,0,30,800,0,10000,sx3event.Energy1,pcevent.Energy2*TMath::Sin(calcsx3theta)); + double sx3theta = TMath::ATan2(sx3rho,sx3z-source_vertex); - 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 pczguess_int = z_to_crossover_rho(pcevent.pos.Z())/TMath::Tan(sx3theta) + source_vertex; + double pcz_guess_int_self = pcevent.pos.Perp()/TMath::Tan(sx3theta) + source_vertex; + plotter->Fill1D("d_guess_self_vs_int_sx3",4000,-200,200,pczguess_int-pcz_guess_int_self); + plotter->Fill2D("guess_self_vs_int_sx3",400,-200,200,400,-200,200,pczguess_int,pcz_guess_int_self); + 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; 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_z_sx3 = x1 + t_minimum*v; + if(vector_closest_to_z_sx3.Perp()>20.0) continue; + double path_length_s = (sx3event.pos-vector_closest_to_z_sx3).Mag()*0.1; + double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1)-path_length_s); + double sinTheta2 = TMath::Sin((sx3event.pos - vector_closest_to_z_sx3).Theta());///TMath::Sin((TVector3(51.5,0,128.) - TVector3(0,0,85)).Theta()); + plotter->Fill2D(std::string("dE3_Ef_Anodesx3B")+"_PC"+std::to_string(phicut),400,0,30,800,0,40000,sx3Efix,pcevent.Energy1*sinTheta2*3.); + plotter->Fill2D(std::string("dE3_E_Anodesx3B")+"_PC"+std::to_string(phicut),400,0,30,800,0,40000,sx3event.Energy1,pcevent.Energy1*sinTheta2*3.); + plotter->Fill2D(std::string("dE3_Ef_Cathodesx3B")+"_PC"+std::to_string(phicut),400,0,30,800,0,10000,sx3Efix,pcevent.Energy2*sinTheta2*3.); + + plotter->Fill1D("VertexReconZ_SX3"+std::to_string(PCSX3TimeCut),600,-1300,1300,vector_closest_to_z_sx3.Z(),"hPCZSX3"); plotter->Fill2D("VertexReconXY_SX3"+std::to_string(PCSX3TimeCut),100,-100,100,100,-100,100,vector_closest_to_z_sx3.X(),vector_closest_to_z_sx3.Y(),"hPCZSX3"); @@ -929,57 +947,74 @@ Bool_t MakeVertex::Process(Long64_t entry) //plotter->Fill2D("pcz_vs_time_strip"+std::to_string(sx3event.ch2),2000,0,2000,600,-200,200,pcevent.Time1*1e-9,pcevent.pos.Z()); //x-axis is all Si det, y-axis is PC anode+cathode only plotter->Fill2D("sx3phi_vs_time",2000,0,2000,180,-360,360,pcevent.Time1*1e-9,sx3event.pos.Phi()*180./M_PI); //x-axis is all Si det, y-axis is PC anode+cathode only + plotter->Fill2D("pcz_vs_sx3pczguess",400,-200,200,400,-200,200,pczguess_int,pcevent.pos.Z()); //x-axis is all Si det, y-axis is PC anode+cathode only + plotter->Fill2D("pcz_vs_sx3pczguess_self",400,-200,200,400,-200,200,pcz_guess_int_self,pcevent.pos.Z()); + plotter->Fill1D("d_pcz_vs_sx3pczguess",400,-200,200,pczguess_int-pcevent.pos.Z()); + plotter->Fill1D("d_pcz_vs_sx3pczguess_self",400,-200,200,pcz_guess_int_self-pcevent.pos.Z()); + + if(pcevent.multi1>0) { //any anodes at all + //plotter->Fill2D("anodeNum_vs_WedgeNum",64,0,64,24,0,24,qqqevent.ch2,i,"onewire"); + } - 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()); + plotter->Fill2D("pcz_vs_sx3pczguess_A1C1",600,-200,200,600,-200,200,pczguess_int,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("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()); - + plotter->Fill2D("pczfix_vs_sx3pczguess_A1C1",600,-200,200,600,-200,200,pczguess_int,pcz_fix_a1c1); + plotter->Fill2D("pcz_vs_sx3pczguess_A1C1_strip"+std::to_string(sx3event.ch2),300,-200,200,600,-200,200,pczguess_int,pcevent.pos.Z()); } - if(pcevent.multi1==1 && pcevent.multi2==2) { - + if(pcevent.multi1==1 && pcevent.multi2==2) { //a1c2 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()); + plotter->Fill2D("pcz_vs_sx3pczguess_A1C2",600,-200,200,600,-200,200,pczguess_int,pcevent.pos.Z()); + plotter->Fill2D("pcz_vs_sx3pczguess_self_A1C2",400,-200,200,400,-200,200,pcz_guess_int_self,pcevent.pos.Z()); + plotter->Fill1D("d_sx3pczguess_minus_pcz_a1c2",400,-200,200,pczguess_int-pcevent.pos.Z()); + plotter->Fill1D("d_sx3pczguess_self_minus_pcz_a1c2",400,-200,200,pcz_guess_int_self-pcevent.pos.Z()); double pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); - + plotter->Fill1D("d_sx3pczguess_minus_pczfix_a1c2",400,-200,200,pczguess_int-pcz_fix); TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix); 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",800,-300,300,r_rhoMin_fix.Z()); + plotter->Fill1D("VertexRecon_pczfix_sx3",800,-300,300,r_rhoMin_fix.Z()); plotter->Fill1D("pczfix_A1C2_1d_sx3",600,-200,200,pcz_fix); - plotter->Fill2D("pczfix_vs_sx3pczguess_A1C2",600,-200,200,600,-200,200,pczguess,pcz_fix); - plotter->Fill2D("pcz_vs_sx3pczguess_A1C2_strip"+std::to_string(sx3event.ch2),300,-200,200,600,-200,200,pczguess,pcevent.pos.Z()); - - double sinTheta_customV = TMath::Sin((sx3event.pos - TVector3(0,0,r_rhoMin_fix.Z())).Theta()); + plotter->Fill2D("pczfix_vs_sx3pczguess_A1C2",600,-200,200,600,-200,200,pczguess_int,pcz_fix); + plotter->Fill2D("pcz_vs_sx3pczguess_A1C2_strip"+std::to_string(sx3event.ch2),300,-200,200,600,-200,200,pczguess_int,pcevent.pos.Z()); + double sinTheta_customV = TMath::Sin((sx3event.pos - r_rhoMin_fix).Theta()); plotter->Fill2D("dE3_E_CathodeSX3_A1C2_TC"+std::to_string(PCSX3TimeCut)+"_PC"+std::to_string(phicut),400,0,30,800,0,10000,sx3event.Energy1,pcevent.Energy2*sinTheta_customV); plotter->Fill2D("dE3_E_AnodeSX3_A1C2_TC"+std::to_string(PCSX3TimeCut)+"_PC"+std::to_string(phicut),400,0,30,800,0,40000,sx3event.Energy1,pcevent.Energy1*sinTheta_customV); - if(TMath::Abs(r_rhoMin_fix.Z())<200.0) { - plotter->Fill2D("dE3_E_AnodeSX3B_A1C2_(vertex_fix_z/100)="+std::to_string(floor(r_rhoMin_fix.Z()/100.0)),400,0,30,800,0,40000,sx3event.Energy1,pcevent.Energy1*sinTheta_customV); - plotter->Fill2D("dE3_E_CathodeSX3B_A1C2_(vertex_fix_z/100)="+std::to_string(floor(r_rhoMin_fix.Z()/100.0)),400,0,30,800,0,10000,sx3event.Energy1,pcevent.Energy2*sinTheta_customV); - } + /*if(TMath::Abs(r_rhoMin_fix.Z())<200.0) { + plotter->Fill2D("dE3_E_AnodeSX3B_A1C2_(vertex_fix_z/100)="+std::to_string(floor(r_rhoMin_fix.Z()/100.0)),400,0,30,800,0,40000,sx3event.Energy1,pcevent.Energy1*sinTheta_customV); + plotter->Fill2D("dE3_E_CathodeSX3B_A1C2_(vertex_fix_z/100)="+std::to_string(floor(r_rhoMin_fix.Z()/100.0)),400,0,30,800,0,10000,sx3event.Energy1,pcevent.Energy2*sinTheta_customV); + }*/ } if(pcevent.multi1==1 && pcevent.multi2==3) { - plotter->Fill2D("pcz_vs_sx3pczguess_A1C3",600,-200,200,600,-200,200,pczguess,pcevent.pos.Z()); - plotter->Fill2D("pcz_vs_sx3pczguess_A1C3_strip"+std::to_string(sx3event.ch2),300,-200,200,600,-200,200,pczguess,pcevent.pos.Z()); + plotter->Fill2D("pcz_vs_sx3pczguess_A1C3",600,-200,200,600,-200,200,pczguess_int,pcevent.pos.Z()); + plotter->Fill2D("pcz_vs_sx3pczguess_A1C3_strip"+std::to_string(sx3event.ch2),300,-200,200,600,-200,200,pczguess_int,pcevent.pos.Z()); } - plotter->Fill2D("pcz_vs_sx3pczguess_int",600,-200,200,600,-200,200,pcz_guess_int,pcevent.pos.Z()); //x-axis is all Si det, y-axis is PC anode+cathode only - plotter->Fill2D("pcz_vs_sx3pczguess_strip"+std::to_string(sx3event.ch2),300,-200,200,600,-200,200,pczguess,pcevent.pos.Z()); + if(pcevent.multi1==2 && pcevent.multi2==1) { + plotter->Fill2D("pcz_vs_sx3pczguess_A2C1",600,-200,200,600,-200,200,pczguess_int,pcevent.pos.Z()); + double pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); + TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix); + 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",800,-300,300,r_rhoMin_fix.Z()); + plotter->Fill1D("pczfix_A2C1_1d_sx3",600,-200,200,pcz_fix); + plotter->Fill2D("pczfix_vs_sx3pczguess_A2C1",600,-200,200,600,-200,200,pczguess_int,pcz_fix); + //plotter->Fill2D("pcz_vs_sx3pczguess_A2C1_strip"+std::to_string(sx3event.ch2),300,-200,200,600,-200,200,pczguess,pcevent.pos.Z()); + } + plotter->Fill2D("pcz_vs_sx3pczguess",600,-200,200,600,-200,200,pczguess_int,pcevent.pos.Z()); //x-axis is all Si det, y-axis is PC anode+cathode only + plotter->Fill2D("pcz_vs_sx3pczguess_strip"+std::to_string(sx3event.ch2),300,-200,200,600,-200,200,pczguess_int,pcevent.pos.Z()); //plotter->Fill2D("pcz_vs_sx3pczguess_phi"+std::to_string(sx3event.pos.Phi()*180/M_PI),300,0,200,600,-200,200,pczguess,pcevent.pos.Z()); - /*plotter->Fill2D("pcz_vs_sx3z_strip="+std::to_string(sx3event.ch2),300,0,100,600,-200,200,sx3z,pcevent.pos.Z(),"sx3_vs_pc_zcorr"); plotter->Fill2D("pcz_vs_sx3z_strip="+std::to_string(sx3event.ch2)+"_a"+std::to_string(pcevent.multi1)+"_c"+std::to_string(pcevent.multi2),300,0,100,600,-200,200,sx3z,pcevent.pos.Z(),"sx3_vs_pc_zcorr"); @@ -1030,49 +1065,58 @@ Bool_t MakeVertex::Process(Long64_t entry) 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 */ - - - for(auto pcevent: PC_Events) { + bool phicut = qqqevent.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/4. && qqqevent.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/4.; 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->Fill2D("dt_pcA_qqqR_vs_qqqRE_PC"+std::to_string(phicut),640,-2000,2000,400,0,30,qqqevent.Time1-pcevent.Time1, qqqevent.Energy1); plotter->Fill1D("dt_pcC_qqqW",640,-2000,2000,qqqevent.Time2 - pcevent.Time2); plotter->Fill2D("phiPC_vs_phiQQQ",180,-360,360,180,-360,360,qqqevent.pos.Phi()*180/M_PI,pcevent.pos.Phi()*180/M_PI); double sinTheta = TMath::Sin((qqqevent.pos - TVector3(0,0,source_vertex)).Theta());///TMath::Sin((TVector3(51.5,0,128.) - TVector3(0,0,85)).Theta()); - + TVector3 x2(pcevent.pos); TVector3 x1(qqqevent.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 r_rhoMin = x1 + t_minimum*v; + if(r_rhoMin.Perp()>20.0) continue; + double path_length_q = (qqqevent.pos-r_rhoMin).Mag()*0.1; + double qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1)-path_length_q); + double sinTheta2 = TMath::Sin((qqqevent.pos - r_rhoMin).Theta());///TMath::Sin((TVector3(51.5,0,128.) - TVector3(0,0,85)).Theta()); + //bool timecut = (qqqevent.Time1 - pcevent.Time1 < -150); - bool timecut = (qqqevent.Time1 - pcevent.Time1 < -150); + bool timecut = (qqqevent.Time1 - pcevent.Time1 < 150); bool lowercut_cath = pcevent.Energy2*sinTheta < 250 && (qqqevent.Energy2 < 5.0 || qqqevent.Energy1 < 5.0) ; - bool phicut = qqqevent.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/4. && qqqevent.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/4.; + if(lowercut_cath && phicut) { plotter->Fill1D("dt_pcA_qqqR_pidlow_PC1",640,-2000,2000,qqqevent.Time1 - pcevent.Time1); plotter->Fill2D("dt_pcA_qqqR_vs_qqqRE_pidlow_PC1",640,-2000,2000,400,0,30,qqqevent.Time1-pcevent.Time1, qqqevent.Energy1); - plotter->Fill1D("dt_pcC_qqqW_pidlow_PC1",640,-2000,2000,qqqevent.Time2 - pcevent.Time2); + plotter->Fill1D("dt_pcC_qqqW_pidlow_PC1",640,-2000,2000,qqqevent.Time2 - pcevent.Time2); } if(timecut) {// && qqqevent.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/4. && qqqevent.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/4. ) { plotter->Fill2D("dE_E_AnodeQQQR",400,0,30,800,0,40000,qqqevent.Energy1,pcevent.Energy1); plotter->Fill2D("dE_E_CathodeQQQR",400,0,30,800,0,10000,qqqevent.Energy2,pcevent.Energy2); + plotter->Fill2D(std::string("dE3_Ef_AnodeQQQR")+"_PC"+std::to_string(phicut),400,0,30,800,0,40000,qqqEfix,pcevent.Energy1*sinTheta2*3.); + plotter->Fill2D(std::string("dE3_Ef_CathodeQQQR")+"_PC"+std::to_string(phicut),400,0,30,800,0,10000,qqqEfix,pcevent.Energy2*sinTheta2*3.); + plotter->Fill2D(std::string("dE3_E_AnodeQQQR")+"_PC"+std::to_string(phicut),400,0,30,800,0,40000,qqqevent.Energy1,pcevent.Energy1*sinTheta2*3.); + plotter->Fill2D(std::string("dE3_E_CathodeQQQR")+"_PC"+std::to_string(phicut),400,0,30,800,0,10000,qqqevent.Energy1,pcevent.Energy2*sinTheta2*3.); + + if(pcevent.multi1==1 && pcevent.multi2==2) plotter->Fill2D("dE_E_AnodeQQQR_a1c2",400,0,30,800,0,40000,qqqevent.Energy1,pcevent.Energy1); if(pcevent.multi1==1 && pcevent.multi2==2) plotter->Fill2D("dE_E_CathodeQQQR_a1c2",400,0,30,800,0,10000,qqqevent.Energy1,pcevent.Energy2); if(pcevent.multi1==2 && pcevent.multi2==1) plotter->Fill2D("dE_E_AnodeQQQR_a2c1",400,0,30,800,0,40000,qqqevent.Energy1,pcevent.Energy1); if(pcevent.multi1==2 && pcevent.multi2==1) plotter->Fill2D("dE_E_CathodeQQQR_a2c1",400,0,30,800,0,10000,qqqevent.Energy1,pcevent.Energy2); - - if(phicut) { + + /*if(phicut) { plotter->Fill2D("dE2_E_AnodeQQQR_TC1PC1_pidlow"+std::to_string(lowercut_cath),400,0,30,800,0,4000,qqqevent.Energy1,pcevent.Energy1*sinTheta); - plotter->Fill2D("dE2_E_CathodeQQQW_TC1PC1_pidlow"+std::to_string(lowercut_cath),400,0,30,800,0,1000,qqqevent.Energy2,pcevent.Energy2*sinTheta); + plotter->Fill2D("dE2_E_CathodeQQQW_TC1PC1_pidlow"+std::to_string(lowercut_cath),400,0,30,800,0,1000,qqqevent.Energy2,pcevent.Energy2*sinTheta); //plotter->Fill2D("E_theta_AnodeQQQR_TC1PC1_pidlow"+std::to_string(lowercut_cath),75,0,90,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1); plotter->Fill2D("E_theta_zoomin_AnodeQQQR_TC1PC1_pidlow"+std::to_string(lowercut_cath),60,0,30,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1); - } - + }*/ + plotter->Fill2D("dE2_E_AnodeQQQR_TC1_PC"+std::to_string(phicut),400,0,30,800,0,4000,qqqevent.Energy1,pcevent.Energy1*sinTheta); plotter->Fill2D("dE2_E_CathodeQQQR_TC1_PC"+std::to_string(phicut),400,0,30,800,0,1000,qqqevent.Energy2,pcevent.Energy2*sinTheta); plotter->Fill2D("dEC_vs_dEA_TC1_PC"+std::to_string(phicut),800,0,40000,800,0,10000,pcevent.Energy1,pcevent.Energy2); @@ -1081,71 +1125,79 @@ Bool_t MakeVertex::Process(Long64_t entry) plotter->Fill1D("dt_pcA_qqqR_timecut",640,-2000,2000,qqqevent.Time1 - pcevent.Time1); plotter->Fill1D("dt_pcC_qqqW_timecut",640,-2000,2000,qqqevent.Time2 - pcevent.Time2); plotter->Fill2D("dE_theta_AnodeQQQR",90,0,90,400,0,20000,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,pcevent.Energy1); - plotter->Fill2D("dE2_theta_AnodeQQQR_zoomin",60,0,30,400,0,5000,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,pcevent.Energy1*sinTheta); + plotter->Fill2D("dE2_theta_AnodeQQQR",90,0,90,400,0,20000,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,pcevent.Energy1*sinTheta); plotter->Fill2D("phiPC_vs_phiQQQ_TimeCut",180,-360,360,180,-360,360,qqqevent.pos.Phi()*180/M_PI,pcevent.pos.Phi()*180/M_PI); - + plotter->Fill1D("d_phiPC_phiQQQ_TimeCut",180,-360,360,qqqevent.pos.Phi()*180/M_PI-pcevent.pos.Phi()*180/M_PI); //plotter->Fill2D("E_theta_AnodeQQQR_TC1_PC"+std::to_string(phicut),75,0,90,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1); //plotter->Fill2D("E_theta_zoomin_AnodeQQQR_TC1_PC"+std::to_string(phicut),60,0,30,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1); //plotter->Fill2D("E2_theta_AnodeQQQR",75,0,90,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1); plotter->Fill2D("Etot2_theta_AnodeQQQR",75,0,90,300,0,15,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,qqqevent.Energy1+pcevent.Energy1*anode_gain*sinTheta); - plotter->Fill2D("dE_theta_CathodeQQQR",75,0,90,800,0,10000,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,pcevent.Energy2); - plotter->Fill2D("dE2_theta_CathodeQQQR",75,0,90,800,0,10000,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,pcevent.Energy2*sinTheta); - plotter->Fill2D("dE2_theta_CathodeQQQR_zoomin",60,0,30,800,0,3000,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,pcevent.Energy2*sinTheta); + //plotter->Fill2D("dE_theta_CathodeQQQR",75,0,90,800,0,10000,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,pcevent.Energy2); + //plotter->Fill2D("dE2_theta_CathodeQQQR",75,0,90,800,0,10000,(qqqevent.pos - TVector3(0,0,source_vertex)).Theta()*180/M_PI,pcevent.Energy2*sinTheta); + - plotter->Fill2D("dE_phi_AnodeQQQR",100,-180,180,800,0,40000,(qqqevent.pos - TVector3(0,0,source_vertex)).Phi()*180/M_PI,pcevent.Energy1); - plotter->Fill2D("dE_phi_CathodeQQQR",100,-180,180,800,0,10000,(qqqevent.pos - TVector3(0,0,source_vertex)).Phi()*180/M_PI,pcevent.Energy2); + //plotter->Fill2D("dE_phi_AnodeQQQR",100,-180,180,800,0,40000,(qqqevent.pos - TVector3(0,0,source_vertex)).Phi()*180/M_PI,pcevent.Energy1); + //plotter->Fill2D("dE_phi_CathodeQQQR",100,-180,180,800,0,10000,(qqqevent.pos - TVector3(0,0,source_vertex)).Phi()*180/M_PI,pcevent.Energy2); plotter->Fill1D("PCZ",800,-200,200,pcevent.pos.Z(),"phicut"); //plotter->Fill1D("PCZ_phicut_a"+std::to_string(aClusters.at(0).size())+"_c"+std::to_string(cClusters.at(0).size()),800,-200,200,pcevent.pos.Z(),"wiremult"); - double pcz_guess_37 = 37./TMath::Tan((qqqevent.pos-TVector3(0,0,source_vertex)).Theta()) + source_vertex; - plotter->Fill2D("pczguess_vs_pc_37",180,0,200,150,0,200,pcz_guess_37,pcevent.pos.Z(),"phicut"); - - double pcz_guess_42 = 42./TMath::Tan((qqqevent.pos-TVector3(0,0,source_vertex)).Theta()) + source_vertex; - plotter->Fill2D("pczguess_vs_pc_42",180,0,200,150,0,200,pcz_guess_42,pcevent.pos.Z(),"phicut"); - double pcz_guess_int = z_to_crossover_rho(pcevent.pos.Z())/TMath::Tan((qqqevent.pos-TVector3(0,0,source_vertex)).Theta()) + source_vertex; + double pcz_guess_int_self = pcevent.pos.Perp()/TMath::Tan((qqqevent.pos-TVector3(0,0,source_vertex)).Theta()) + source_vertex; + plotter->Fill1D("d_guess_self_vs_int_qqq",400,-200,200,pcz_guess_int-pcz_guess_int_self); + plotter->Fill2D("guess_self_vs_int_qqq",400,-200,200,400,-200,200,pcz_guess_int,pcz_guess_int_self); + + + //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"); + plotter->Fill2D("pczguess_vs_pc_int_self",400,-200,200,600,-200,200,pcz_guess_int_self,pcevent.pos.Z(),"phicut"); + plotter->Fill1D("d_pczqqq_vs_pc_int",400,-200,200,pcz_guess_int-pcevent.pos.Z()); + plotter->Fill1D("d_pczqqq_vs_pc_int_self",400,-200,200,pcz_guess_int_self-pcevent.pos.Z()); 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); + 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) { + if(pcevent.multi1==1 && pcevent.multi2==2) { //a1c2 double pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix); 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; - double sinTheta_customV = TMath::Sin((qqqevent.pos - TVector3(0,0,r_rhoMin_fix.Z())).Theta()); + double sinTheta_customV = TMath::Sin((qqqevent.pos - r_rhoMin_fix).Theta()); plotter->Fill2D("dE3_E_CathodeQQQW_A1C2_TC1_PC"+std::to_string(phicut),400,0,30,800,0,10000,qqqevent.Energy2,pcevent.Energy2*sinTheta_customV); plotter->Fill2D("dE3_E_AnodeQQQR_A1C2_TC1_PC"+std::to_string(phicut),400,0,30,800,0,10000,qqqevent.Energy1,pcevent.Energy1*sinTheta_customV); - + plotter->Fill1D("VertexRecon_pczfix_qqq",800,-400,400,r_rhoMin_fix.Z()); plotter->Fill1D("VertexRecon_pczfix_qqq_PC"+std::to_string(phicut)+"_pidlow"+std::to_string(lowercut_cath),800,-400,400,r_rhoMin_fix.Z()); - - if(TMath::Abs(r_rhoMin_fix.Z())<200.0) { + plotter->Fill1D("d_qqqpczguess_minus_pcz_a1c2",400,-200,200,pcz_guess_int-pcevent.pos.Z()); + plotter->Fill1D("d_qqqpczguess_self_minus_pcz_a1c2",400,-200,200,pcz_guess_int_self-pcevent.pos.Z()); + /*if(TMath::Abs(r_rhoMin_fix.Z())<200.0) { plotter->Fill2D("dE3_E_AnodeQQQR_A1C2_(vertex_fix_z/100)="+std::to_string(floor(r_rhoMin_fix.Z()/100.0)),400,0,30,800,0,40000,qqqevent.Energy1,pcevent.Energy1*sinTheta_customV); plotter->Fill2D("dE3_E_CathodeQQQR_A1C2_(vertex_fix_z/100)="+std::to_string(floor(r_rhoMin_fix.Z()/100.0)),400,0,30,800,0,10000,qqqevent.Energy1,pcevent.Energy2*sinTheta_customV); - } + } */ - plotter->Fill1D("pczfix_A1C2_1d_qqq",600,-200,200,pcz_fix); - plotter->Fill2D("pczfix_vs_qqqpczguess_A1C2",600,-200,200,600,-200,200,pcz_guess_int,pcz_fix); + plotter->Fill1D("pczfix_A1C2_1d_qqq",600,-200,200,pcz_fix); + plotter->Fill1D("d_qqqpczguess_minus_pczfix_a1c2",400,-200,200,pcz_guess_int-pcz_fix); + plotter->Fill1D("d_qqqpczguess_self_minus_pczfix_a1c2",400,-200,200,pcz_guess_int_self-pcz_fix); + plotter->Fill2D("pczfix_vs_qqqpczguess_A1C2",600,-200,200,600,-200,200,pcz_guess_int,pcz_fix); + plotter->Fill2D("pczfix_vs_qqqpczguess_self_A1C2",600,-200,200,600,-200,200,pcz_guess_int_self,pcz_fix); plotter->Fill2D("pczguess_vs_pc_int_A1C2",400,-200,200,600,-400,400,pcz_guess_int,pcevent.pos.Z(),"phicut"); - + plotter->Fill2D("pczguess_self_vs_pc_int_A1C2",400,-200,200,600,-400,400,pcz_guess_int_self,pcevent.pos.Z(),"phicut"); + double path_length = (qqqevent.pos-TVector3(0,0,r_rhoMin_fix.Z())).Mag()*0.1; //std::cout << path_length << std::endl; double qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1)-path_length); double qqqEfix_p = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1)-path_length); - + plotter->Fill2D("E_thetaf_AnodeQQQR_TC1_PC"+std::to_string(phicut),180,0,180,600,0,15,(qqqevent.pos - TVector3(0,0,r_rhoMin_fix.Z())).Theta()*180/M_PI,qqqevent.Energy1); if(lowercut_cath) plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC"+std::to_string(phicut)+"_pidlow"+std::to_string(lowercut_cath),180,0,180,600,0,15,(qqqevent.pos - TVector3(0,0,r_rhoMin_fix.Z())).Theta()*180/M_PI,qqqEfix_p); @@ -1153,449 +1205,303 @@ Bool_t MakeVertex::Process(Long64_t entry) std::string zcut = "_"+std::to_string((TMath::Abs(r_rhoMin_fix.Z())<180)); plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC"+std::to_string(phicut)+"_pidlow"+std::to_string(lowercut_cath)+zcut,180,0,180,600,0,15,(qqqevent.pos - TVector3(0,0,r_rhoMin_fix.Z())).Theta()*180/M_PI,qqqEfix); } - std::string morecuts = "_pidlow"+std::to_string(lowercut_cath)+"_vertexfix="+std::to_string(floor(r_rhoMin_fix.Z()/20)*20+10); //plotter->Fill2D("E_thetaf_AnodeQQQR_TC1_PC"+std::to_string(phicut)+morecuts,180,0,180,800,0,8,(qqqevent.pos - TVector3(0,0,r_rhoMin_fix.Z())).Theta()*180/M_PI,qqqevent.Energy1,"morecuts"); - //plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC"+std::to_string(phicut)+morecuts,180,0,180,800,0,8,(qqqevent.pos - TVector3(0,0,r_rhoMin_fix.Z())).Theta()*180/M_PI,qqqEfix,"morecuts"); - plotter->Fill2D("dE3_Ef_AnodeQQQR_TC1"+std::to_string(phicut)+"_pidlow"+std::to_string(lowercut_cath),600,0,15,800,0,40000,qqqEfix,pcevent.Energy1*sinTheta_customV); plotter->Fill2D("dE3_Ef_CathodeQQQR_TC1PC"+std::to_string(phicut)+"_pidlow"+std::to_string(lowercut_cath),600,0,15,800,0,10000,qqqEfix,pcevent.Energy2*sinTheta_customV); } + if(pcevent.multi1==2 && pcevent.multi2==1) { + plotter->Fill2D("pcz_int_vs_qqqpczguess_A2C1",600,-200,200,600,-200,200,pcz_guess_int,pcevent.pos.Z()); + plotter->Fill2D("pczguess_vs_pc_int_self_A2C1",400,-200,200,600,-200,200,pcz_guess_int_self,pcevent.pos.Z(),"phicut"); + double pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); + TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix); + 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_a2c1_pczfix_qqq",800,-300,300,r_rhoMin_fix.Z()); + plotter->Fill1D("pczfix_A2C1_1d_qqqs",600,-200,200,pcz_fix); + plotter->Fill2D("pczfix_vs_qqqpczguess_A2C1",600,-200,200,600,-200,200,pcz_guess_int,pcz_fix); + //plotter->Fill2D("pcz_vs_sx3pczguess_A2C1_strip"+std::to_string(sx3event.ch2),300,-200,200,600,-200,200,pczguess,pcevent.pos.Z()); + } + double qqqrho = qqqevent.pos.Perp(); double qqqz = (qqqevent.pos - TVector3(0,0,source_vertex)).Z(); double tan_theta = qqqrho/qqqz; double pcz_guess_int2 = z_to_crossover_rho(pcevent.pos.Z())/tan_theta + source_vertex; plotter->Fill2D("pczguess_vs_pc_int2",180,0,200,150,0,200,pcz_guess_int2,pcevent.pos.Z(),"phicut"); - + double qqqz2 = (qqqevent.pos - r_rhoMin).Z(); double tan_theta2 = qqqrho/qqqz2; double pcz_guess_int3 = z_to_crossover_rho(pcevent.pos.Z())/tan_theta2 + r_rhoMin.Z(); plotter->Fill2D("pczguess_vs_pc_int3",180,0,200,150,0,200,pcz_guess_int3,pcevent.pos.Z(),"phicut"); //plotter->Fill2D("pczguess_vs_pc_int2_a"+std::to_string(pcevent.multi1)+"_c"+std::to_string(pcevent.multi2),180,0,200,150,0,200,pcz_guess_int2,pcevent.pos.Z(),"phicut"); - double pcz_guess = pcz_guess_int; plotter->Fill2D("pctheta_vs_qqqtheta_sv",180,-360,360,180,-360,360,(qqqevent.pos-TVector3(0,0,source_vertex)).Theta()*180/M_PI,(pcevent.pos-TVector3(0,0,source_vertex)).Theta()*180/M_PI,"phicut"); plotter->Fill2D("pctheta_vs_qqqtheta_rmz",180,-360,360,180,-360,360,(qqqevent.pos-TVector3(0,0,r_rhoMin.Z())).Theta()*180/M_PI,(pcevent.pos-TVector3(0,0,r_rhoMin.Z())).Theta()*180/M_PI,"phicut"); plotter->Fill2D("pctheta_vs_qqqtheta_rm",180,-360,360,180,-360,360,(qqqevent.pos-r_rhoMin).Theta()*180/M_PI,(pcevent.pos-r_rhoMin).Theta()*180/M_PI,"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->Fill2D("pczguess_vs_pc_phi="+std::to_string(qqqevent.pos.Phi()*180./M_PI),300,0,200,150,0,200,pcz_guess_int,pcevent.pos.Z(),"phicut"); } } }//end PC QQQ coincidence - //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_Cathode_Coincidence_Matrix_qqq"+std::to_string(HitNonZero), 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(cID, cE)); - cESum += cE; - } - }*/ - if((aIDMax + cID)%24 == 22 || (aIDMax + cID)%24==23 || (aIDMax + cID)%24>=0 || (aIDMax + cID)%24<=3 ) { - corrcatMax.push_back(std::pair(cID, cE)); - cESum += cE; - if(cE > cEMax) { - cEMax = cE; - cIDMax = cID; - } - } - } - } - } - - 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 (pwinstance.Crossover[aIDMax][corr.first][0].z > 9000000) - continue; - if (cESum > 0) - { - x += (corr.second) / cESum * pwinstance.Crossover[aIDMax][corr.first][0].x; - y += (corr.second) / cESum * pwinstance.Crossover[aIDMax][corr.first][0].y; - z += (corr.second) / cESum * pwinstance.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); - if(realtime) { - //crossoverg->SetPoint(0,x,y,z); - crossoverg->AddPoint(x,y,z); - } - //std::cout << "Anode Intersection: " << anodeIntersection.X() << ", " << anodeIntersection.Y() << ", " << anodeIntersection.Z() << " " << aIDMax << 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; - } - - if(anodeIndex!=-1 && cathodeIndex !=-1 && hitPos.Perp()!=0 && anodeIntersection.Perp()!=0 && realtime && PCQQQPhiCut && PCQQQTimeCut) { - //can1->Modified(); - //can1->Update(); - TVector3 x2(anodeIntersection); - TVector3 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()); - TVector3 r_rhoMin = x1 + t_minimum*v; - - trajectory->SetPoint(0,x1.X(),x1.Y(),x1.Z()); - trajectory->SetPoint(1,r_rhoMin.X(),r_rhoMin.Y(),r_rhoMin.Z()); - - for(auto cath: corrcatMax) { - plc[cath.first]->SetLineWidth(3); - //plc[cath.first]->SetLineStyle(kLine); - } - for(auto anodeW: anodeHits) { - pla[anodeW.first]->SetLineWidth(3); - //pla[anodeW.first]->SetLineStyle(kLine); - } - //can2->Modified(); - //can2->Update(); - //while(can1->WaitPrimitive()); - - //pla[anodeIndex]->SetLineWidth(1); - //pla[anodeIndex]->SetLineStyle(kDotted); - for(auto anodeW: anodeHits) { - pla[anodeW.first]->SetLineWidth(1); - pla[anodeW.first]->SetLineStyle(kDotted); - } - for(auto cath: corrcatMax) { - plc[cathodeIndex]->SetLineStyle(kDotted); - plc[cath.first]->SetLineWidth(1); - } - } - - 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 && anodeIntersection.Perp() > 0 && PCSX3TimeCut) { - plotter->Fill1D("PC_Z_Projection_sx3", 600, -200, 200, anodeIntersection.Z(), "hPCZSX3"); - } - 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(qqqenergy<4.0) - plotter->Fill1D("VertexRecon_Z(qqqE<4.0MeV)_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,20, 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_TimeCut" + std::to_string(qqq.id[i]), 400, -100, 100, 400, -100, 100,hitPos.X(),hitPos.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(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(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"); - plotter->Fill2D("PC_Z_vs_QQQRho", 600, -300, 300, 40, 40, 110, anodeIntersection.Z(), hitPos.Perp(), "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(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,source_vertex); //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 = 37.0/TMath::Tan((hitPos-guessVertex).Theta()) + guessVertex.Z(); //this is ideally kept to be all QQQ+userinput for calibration of pcz - double pcz_guess = z_to_crossover_rho(anodeIntersection.Z())/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=0.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),200,0,200,200,0,200,pcz_guess,anodeIntersection.Z()+pczoffset,"pczguess"); //entirely qqq-derived position vs entirely PC derived position - plotter->Fill2D("pczguess_vs_pcz",200,0,200,200,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"); - } - - if(anodeIntersection.Perp()!=0) { - plotter->Fill2D("AnodeMaxE_Vs_Cathode_Sum_Energy", 2000, 0, 20000, 2000, 0, 10000, aEMax, cESum, "hGMPC"); - plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy", 800, 0, 20000, 800, 0, 10000, aESum, cEMax, "hGMPC"); - plotter->Fill2D("AnodeMaxE_Vs_Cathode_Max_Energy", 800, 0, 20000, 800, 0, 10000, aEMax, cEMax, "hGMPC"); - //double sinTheta = TMath::Sin((anodeIntersection - TVector3(0,0,source_vertex)).Theta());///TMath::Sin((TVector3(51.5,0,128.) - TVector3(0,0,85)).Theta()); - //plotter->Fill2D("AnodeMaxE_Vs_Cathode_Max_Energy_path_corrected", 800, 0, 20000, 800, 0, 10000, aEMax*sinTheta, cEMax*sinTheta, "hGMPC"); - plotter->Fill2D("AnodeSumE_Vs_Cathode_Sum_Energy", 800, 0, 20000, 800, 0, 10000, aESum, cESum, "hGMPC"); - plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_TC"+std::to_string(PCQQQTimeCut)+"_PC"+std::to_string(PCQQQPhiCut), 800, 0, 20000, 800, 0, 10000, aESum, cEMax, "hGMPC"); - //plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_path_corrected"+std::to_string(PCQQQTimeCut)+"_PC"+std::to_string(PCQQQPhiCut), 800, 0, 20000, 800, 0, 10000, aESum*sinTheta, cEMax*sinTheta, "hGMPC"); - //plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_path_corrected", 800, 0, 20000, 800, 0, 10000, aESum*sinTheta, cEMax*sinTheta, "hGMPC"); - - if(PCQQQTimeCut && PCQQQPhiCut) { - plotter->Fill2D("AnodeSumE_Vs_Cathode_Max_Energy_TC"+std::to_string(PCQQQTimeCut)+"_PC"+std::to_string(PCQQQPhiCut)+"_cMax"+std::to_string(cIDMax), 800, 0, 20000, 800, 0, 10000, aESum, cEMax, "hGMPC"); - } - //plotter->Fill2D("AnodeSumE_Vs_CathodeSum_Energy_path_corrected", 800, 0, 20000, 800, 0, 10000, aESum*sinTheta, cESum*sinTheta, "hGMPC"); - //plotter->Fill2D("AnodeSumE_Vs_CathodeSum_Energy_path_corrected_TC"+std::to_string(PCQQQTimeCut)+"_PC"+std::to_string(PCQQQPhiCut), 800, 0, 20000, 800, 0, 10000, aESum*sinTheta, cESum*sinTheta, "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"); - } - - for(auto cwevent: cWireEvents) { - //plotter->Fill1D("cwdtqqq_vs_cw"+std::to_string(PCQQQTimeCut),800,-2000,2000,24,0,24,std::get<2>(cwevent)-qqqtimestamp,std::get<0>(cwevent)); - for(auto awevent: aWireEvents) { - plotter->Fill2D("aw_vs_cw",24,0,24,24,0,24,std::get<0>(awevent),std::get<0>(cwevent)); - plotter->Fill2D("aw_vs_cw_dtq"+std::to_string(PCQQQTimeCut),24,0,24,24,0,24,std::get<0>(awevent),std::get<0>(cwevent)); - } - } - for(auto awevent: aWireEvents) { - //plotter->Fill1D("awdtqqq_vs_aw"+std::to_string(PCQQQTimeCut),800,-2000,2000,24,0,24,std::get<2>(awevent)-qqqtimestamp,std::get<0>(awevent)); - } - return kTRUE; + return kTRUE; } void MakeVertex::Terminate() { plotter->FlushToDisk(10); -/* can1->Modified(); - can1->Update(); - can2->Modified(); - can2->Update(); - while(can1->WaitPrimitive()); - while(can2->WaitPrimitive());*/ +} + +void miscHistograms_oneWire(HistPlotter* plotter, std::vector QQQ_Events, std::vector>> aClusters) { + //consider the 'proton-like' QQQ branch seen in a,p data + TRandom3 rand; + rand.SetSeed();//random seed set + Kinematics apkin_a(1.008664916,4.002603254,4.002603254,1.008664916,7.1); //m3 is alpha, 6.79 MeV is 7.0 MeV proton energy after kapton+100mm 4He gas (molar mass 5.2, 250 torr) + for(auto qqqevent: QQQ_Events) { + if(qqqevent.Energy1 < 0.6) continue; //coarse gating + //if(qqqevent.Energy1 > 5.0) continue; //coarse gating + for(const auto acluster: aClusters) { + auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(acluster,"ANODE"); + //if(apSumE<6000) continue; + int a_number = acluster.size(); + TVector3 pc_closest = pwinstance.getClosestWirePosAtWirePhi(apwire,qqqevent.pos.Phi()); + plotter->Fill1D("dt_anode_interp_qqq",800,-2000,2000,qqqevent.Time1 - apTSMaxE,"ainterp_noc"); + if(qqqevent.Time1 - apTSMaxE < 150) { + bool phicut = qqqevent.pos.Phi() <= pc_closest.Phi()+TMath::Pi()/4. && qqqevent.pos.Phi() >= pc_closest.Phi()-TMath::Pi()/4.; + + TVector3 x2(pc_closest), x1(qqqevent.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 r_rhoMin_fix = x1 + t_minimum*v; + + double theta_q = (qqqevent.pos - r_rhoMin_fix).Theta(); + double sinTheta2 = TMath::Sin(theta_q); + + if(r_rhoMin_fix.Perp()>6.0) continue; + if(r_rhoMin_fix.Z() < -173.6 || r_rhoMin_fix.Z() > 100) continue; + if(!phicut) continue; + plotter->Fill1D("dt_anode_ainterp_qqq_gated",800,-2000,2000,qqqevent.Time1 - apTSMaxE,"ainterp_noc"); + plotter->Fill2D("dt_anode_ainterp_qqq_gated_vs_qqqE",800,-2000,2000,800,0,10,qqqevent.Time1 - apTSMaxE,qqqevent.Energy1,"ainterp_noc"); + plotter->Fill2D("dEa_ainterp_Eqqq_TC1_ignC_a"+std::to_string(acluster.size()),400,0,10,800,0,40000,qqqevent.Energy1,apSumE,"ainterp_noc"); + plotter->Fill2D("pcPhi_ainterp_qqqPhi_TC1_ignC_a"+std::to_string(acluster.size()),120,-360,360,120,-360,360,pc_closest.Phi()*180./M_PI,qqqevent.pos.Phi()*180./M_PI,"ainterp_noc"); + plotter->Fill2D("pcZ_ainterp_qqqZ_TC1_ignC_a"+std::to_string(acluster.size())+"_PC"+std::to_string(phicut),300,-100,200,400,-200,200,qqqevent.pos.Z(),pc_closest.Z(),"ainterp_noc"); + + + //plotter->Fill2D("pcZ_ainterp_qqqpczguess_TC1_ignC_a"+std::to_string(acluster.size()),300,-100,200,400,-200,200,pczguess,pc_closest.Z(),"ainterp_noc"); + plotter->Fill2D("dEa3_ainterp_Eqqq_TC1_ignC_a"+std::to_string(acluster.size())+"_PC"+std::to_string(phicut),1200,0,30,800,0,30000,qqqevent.Energy1,apSumE*sinTheta2*3.,"ainterp_noc"); + + plotter->Fill2D("vertexZ_ainterp_qqqZ_TC1_ignC_a"+std::to_string(acluster.size()),300,-100,200,800,-400,400,qqqevent.pos.Z(),r_rhoMin_fix.Z(),"ainterp_noc"); + plotter->Fill1D("vertexZ1d_ainterp_qqqZ_TC1_ignC_a"+std::to_string(acluster.size()),800,-400,400,r_rhoMin_fix.Z(),"ainterp_noc"); + plotter->Fill2D("vertexXY_ainterp_TC1_ignC_a"+std::to_string(acluster.size()),200,-100,100,200,-100,100,r_rhoMin_fix.X(),r_rhoMin_fix.Y(),"ainterp_noc"); + + double path_length_q = (qqqevent.pos-r_rhoMin_fix).Mag()*0.1; + double qqqEfix; + qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1)-path_length_q); + plotter->Fill1D("pmisc_ow_Ex_from_alpha",200,-10,10,apkin_a.getExc(qqqEfix,theta_q*180/M_PI),"ainterp_noc"); + plotter->Fill2D("pmisc_ow_Ef_vs_theta_qqq",100,0,180,800,0,20,theta_q*180/M_PI,qqqEfix,"ainterp_noc"); + plotter->Fill2D("pmisc_ow_VertexReconZ_vs_Ef",800,-400,400,800,0,20,r_rhoMin_fix.Z(),qqqEfix,"ainterp_noc"); + } + } + /* + for(int i=0; i<24; i++) { + if(aWireEvents.find(i) != aWireEvents.end()) { + auto awire = aWireEvents[i]; + if(qqqevent.Time1 -(double)std::get<2>(awire)< 150) { + //plotter->Fill2D("qqq_z_phi2_awire"+std::to_string(std::get<0>(awire)), 400,-100,100, 100, -200,200,qqqevent.pos.Z(), qqqevent.pos.Phi()*180/M_PI ); + //plotter->Fill2D("qqq_z_strip#_awire"+std::to_string(std::get<0>(awire)), 400,-100,100, 100, -50,50,qqqevent.pos.Z(), qqqevent.ch2); + plotter->Fill2D("anodeNum_vs_stripNum",64,0,64,24,0,24,qqqevent.ch2,i,"onewire"); + bool qqqdiagonalphi = (!plotter->FindCut("anode_qqq_diag1")->IsInside(qqqevent.ch2,i)) || (!plotter->FindCut("anode_qqq_diag2")->IsInside(qqqevent.ch2,i)); + plotter->Fill2D("anodeNum_vs_stripNum_diag"+std::to_string(qqqdiagonalphi),64,0,64,24,0,24,qqqevent.ch2,i,"onewire"); + plotter->Fill2D("onewire_dEa_Eqqq_TC1_fullev"+std::to_string(PC_Events.size()>0)+"_PC"+std::to_string(qqqdiagonalphi),400,0,10,800,0,40000,qqqevent.Energy1,std::get<1>(awire),"onewire"); + //plotter->Fill2D("onewire_aNum_qqqPhi_TC1_fullev"+std::to_string(PC_Events.size()>0),24,0,24,120,-360,360,i,qqqevent.pos.Phi()*180./M_PI,"onewire"); + //plotter->Fill2D("qqq_z_phi_ow_awire"+std::to_string(anodeIndex)+"_qqqstrip"+std::to_string(qqqevent.ch2), 400,-100,100, 200, -200,200,qqqevent.pos.Z(), qqqevent.pos.Phi()*180/M_PI ); + } + } + + if(cWireEvents.find(i) != cWireEvents.end()) { + auto cwire = cWireEvents[i]; + if(qqqevent.Time1 -(double)std::get<2>(cwire) < 150) { + //plotter->Fill2D("qqq_z_phi2_cwire"+std::to_string(std::get<0>(cwire)),400,-100,100, 100, -200,200,qqqevent.pos.Z(), qqqevent.pos.Phi()*180/M_PI ); + //plotter->Fill2D("qqq_z_strip#_cwire"+std::to_string(std::get<0>(cwire)),400,-100,100, 100, -50,50,qqqevent.pos.Z(), qqqevent.ch2 ); + plotter->Fill2D("onewire_dEc_Eqqq_fullev"+std::to_string(PC_Events.size()>0),400,0,10,800,0,40000,qqqevent.Energy1,std::get<1>(cwire),"onewire"); + plotter->Fill2D("onewire_cNum_qqqPhi_TC1_fullev"+std::to_string(PC_Events.size()>0),24,0,24,120,-360,360,i,qqqevent.pos.Phi()*180./M_PI,"onewire"); + } + } + }//for 'i' loop*/ + }//end QQQEvents loop +} + + +void protonMiscHistograms(HistPlotter* plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events) { + //consider the 'proton-like' QQQ branch seen in a,p data + TRandom3 rand; + rand.SetSeed();//random seed set + Kinematics apkin_a(1.008664916,4.002603254,4.002603254,1.008664916,7.1); //m3 is alpha, 6.79 MeV is 7.0 MeV proton energy after kapton+100mm 4He gas (molar mass 5.2, 250 torr) + for(auto qqqevent: QQQ_Events) { + if(qqqevent.Energy1 < 0.6) continue; //coarse gating + //if(qqqevent.Energy1 > 5.0) continue; //coarse gating + for(auto pcevent: PC_Events) { + if(!(pcevent.multi1==1 && pcevent.multi2<=2)) continue; + //if(pcevent.Energy1 > 11000) continue; //coarse gating + bool phicut = qqqevent.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/4. && qqqevent.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/4.; + if(!phicut) continue; + //if(pcevent.Time1-qqqevent.Time1<-150 || pcevent.Time1-qqqevent.Time1 >850) continue; + + double pcz_fix, pcz_dith=pcevent.pos.Z(); + if(pcevent.multi2==2) + pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); + else { + pcz_fix = rand.Gaus(pcevent.pos.Z(),8.0);//dither for a1c1 events + pcz_dith = pcz_fix; + } + TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix); + TVector3 x1(qqqevent.pos); + 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; + double vertex_z = r_rhoMin_fix.Z(); + //double theta_q = (qqqevent.pos - TVector3(0,0,vertex_z)).Theta(); + double theta_q = (qqqevent.pos - r_rhoMin_fix).Theta(); + double sinTheta_customV = TMath::Sin(theta_q); + //if(r_rhoMin_fix.Perp()>6) continue; + bool cathode_alpha_select = (pcevent.Energy2 > 1400); + if(vertex_z < -173.6 || vertex_z > 100) continue; + + + //What's below: radial cut, time coincident, phi-correlated events with possible energy selection applied to both E_si and dE_Anodes + auto plot_with_tag = [&](std::string tag="") { + std::string pmlabel = "proton+misc"+tag; + plotter->Fill2D("pmisc_dE_E_AnodeQQQ"+tag,400,0,10,800,0,40000,qqqevent.Energy1,pcevent.Energy1,pmlabel); + plotter->Fill2D("pmisc_dE_E_CathodeQQQ"+tag,400,0,10,800,0,10000,qqqevent.Energy1,pcevent.Energy2,pmlabel); + plotter->Fill2D("pmisc_dE3_E_AnodeQQQ"+tag,400,0,10,400,0,40000,qqqevent.Energy1,pcevent.Energy1*sinTheta_customV*3.,pmlabel); + plotter->Fill2D("pmisc_dE3_E_CathodeQQQ"+tag,400,0,10,400,0,10000,qqqevent.Energy1,pcevent.Energy2*sinTheta_customV,pmlabel); + plotter->Fill2D("pmisc_dPhi_QQQ_PC"+tag,180,-360,360,180,-360,360,pcevent.pos.Phi()*180/M_PI,qqqevent.pos.Phi()*180/M_PI,pmlabel); + plotter->Fill1D("pmisc_dt_Anode_QQQ_PC"+std::to_string(phicut)+tag,600,-2000,2000,pcevent.Time1-qqqevent.Time1,pmlabel); + plotter->Fill1D("pmisc_dt_Cathode_QQQ"+tag,600,-2000,2000,pcevent.Time2-qqqevent.Time1,pmlabel); + plotter->Fill2D("pmisc_dt_Anode_E_QQQ_PC"+std::to_string(phicut)+tag,600,-2000,2000,400,0,10,pcevent.Time1-qqqevent.Time1,qqqevent.Energy1,pmlabel); + plotter->Fill2D("pmisc_dt_AnodeQQQ_vsPCPhi"+tag,600,-2000,2000,180,-360,360,pcevent.Time1-qqqevent.Time1,pcevent.pos.Phi()*180./M_PI,pmlabel); + plotter->Fill2D("pmisc_dt_Cathode_E_QQQ"+tag,600,-2000,2000,400,0,10,pcevent.Time2-qqqevent.Time1,qqqevent.Energy1,pmlabel); + plotter->Fill2D("pmisc_dt_CathodeQQQ_vsPCPhi"+tag,600,-2000,2000,180,-360,360,pcevent.Time2-qqqevent.Time1,pcevent.pos.Phi()*180./M_PI,pmlabel); + plotter->Fill1D("pmisc_pczfix"+tag,600,-300,300,pcz_fix,pmlabel); + if(pcevent.multi2==2) { + plotter->Fill1D("pmisc_pcz"+tag,600,-300,300,pcevent.pos.Z(),pmlabel); + plotter->Fill1D("pmisc_pcz2"+tag,600,-300,300,pcevent.pos.Z(),pmlabel); + } + if(pcevent.multi2==1) { + plotter->Fill1D("pmisc_pcz"+tag,600,-300,300,pcz_fix,pmlabel); + plotter->Fill1D("pmisc_pcz1"+tag,600,-300,300,pcevent.pos.Z(),pmlabel); + plotter->Fill1D("pmisc_pcz_dith"+tag,600,-300,300,pcz_dith,pmlabel); + } + + //double path_length_q = (qqqevent.pos-TVector3(0,0,vertex_z)).Mag()*0.1; + //double path_length_s = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; + double path_length_q = (qqqevent.pos-r_rhoMin_fix).Mag()*0.1; + double qqqEfix; + if(tag == "_cathode_alphas") {//satisfied when find succeeds + qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1)-path_length_q); + plotter->Fill1D("pmisc_Ex_from_alpha",200,-10,10,apkin_a.getExc(qqqEfix,theta_q*180/M_PI),pmlabel); + } + else + qqqEfix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(qqqevent.Energy1)-path_length_q); + //plotter->Fill2D("qqqEf_sx3E_matrix_all"+tag,400,0,10,400,0,10,qqqEfix,sx3event.Energy1,pmlabel); + plotter->Fill2D("pmisc_dE3_Ef_AnodeQQQ"+tag,400,0,10,400,0,40000,qqqEfix,pcevent.Energy1*sinTheta_customV*3,pmlabel); + plotter->Fill2D("pmisc_dE3_Ef_CathodeQQQ"+tag,400,0,10,400,0,10000,qqqEfix,pcevent.Energy2*sinTheta_customV,pmlabel); + + plotter->Fill1D("pmisc_VertexReconZ"+tag,800,-400,400,vertex_z,pmlabel); + plotter->Fill2D("pmisc_VertexReconXY"+tag,200,-100,100,200,-100,100,r_rhoMin_fix.X(),r_rhoMin_fix.Y(),pmlabel); + plotter->Fill2D("pmisc_VertexReconZ_vs_Ef"+tag,800,-400,400,800,0,20,vertex_z,qqqEfix,pmlabel); + plotter->Fill2D("pmisc_VertexReconZ_vs_Ef"+tag+"_a"+std::to_string(pcevent.multi1),800,-400,400,800,0,20,vertex_z,qqqEfix,pmlabel); + + plotter->Fill2D("pmisc_Ef_vs_theta_qqq"+tag,100,0,180,800,0,20,theta_q*180/M_PI,qqqEfix,pmlabel); + if(pcevent.multi2==1) { + plotter->Fill2D("pmisc_Ef_vs_theta_qqq_a1c1"+tag,100,0,180,800,0,20,theta_q*180/M_PI,qqqEfix,pmlabel); + plotter->Fill2D("pmisc_VertexReconZ_vs_Ef_a1c1"+tag,800,-400,400,800,0,20,vertex_z,qqqEfix,pmlabel); + } + }; + + if(cathode_alpha_select) + plot_with_tag("_cathode_alphas"); + //else + // plot_with_tag("_cathode_protons"); + //plot_with_tag(); + + //plotter->Fill1D("pmisc_Ex_from_protons",200,-10,10,apkin_p.getExc(qqqEfix,theta_s*180/M_PI),pmlabel); + + }//end PCEvents loop + }//end QQQEvents loop +} + +void protonMiscHistograms_sx3(HistPlotter* plotter, std::vector QQQ_Events, std::vector SX3_Events, std::vector PC_Events) { + //consider the 'proton-like' QQQ branch seen in a,p data + for(auto sx3event: SX3_Events) { + if(sx3event.Energy1 < 1.2) continue; //coarse gating + //if(sx3event.Energy1 > 5.0) continue; //coarse gating + for(auto pcevent: PC_Events) { + if(!(pcevent.multi1==1 && pcevent.multi2==2)) continue; + //if(pcevent.Energy1 > 11000) continue; //coarse gating + bool phicut = sx3event.pos.Phi() <= pcevent.pos.Phi()+TMath::Pi()/3. && sx3event.pos.Phi() >= pcevent.pos.Phi()-TMath::Pi()/3.; + if(!phicut) continue; + //if(pcevent.Time1-sx3event.Time1<-150 || pcevent.Time1-sx3event.Time1 >850) continue; + + double pcz_fix = pcfix_func.Eval(pcevent.pos.Z()); + TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix); + TVector3 x1(sx3event.pos); + 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; + double vertex_z = r_rhoMin_fix.Z(); + //double theta_q = (sx3event.pos - TVector3(0,0,vertex_z)).Theta(); + + if(r_rhoMin_fix.Perp()>10.0) continue; + double theta_s = (sx3event.pos - r_rhoMin_fix).Theta(); + double sinTheta_customV = TMath::Sin(theta_s); + bool cathode_alpha_select = (pcevent.Energy2 > 1400); + + + + //What's below: radial cut, time coincident, phi-correlated events with possible energy selection applied to both E_si and dE_Anodes + auto plot_with_tag = [&](std::string tag="") { + std::string pmlabel = "proton+miscsx3"+tag; + plotter->Fill2D("pmiscs_dE_E_Anodesx3"+tag,400,0,10,800,0,40000,sx3event.Energy1,pcevent.Energy1,pmlabel); + plotter->Fill2D("pmiscs_dE_E_Cathodesx3"+tag,400,0,10,800,0,10000,sx3event.Energy1,pcevent.Energy2,pmlabel); + plotter->Fill2D("pmiscs_dE3_E_Anodesx3"+tag,400,0,10,400,0,40000,sx3event.Energy1,pcevent.Energy1*sinTheta_customV*3.,pmlabel); + plotter->Fill2D("pmiscs_dE3_E_Cathodesx3"+tag,400,0,10,400,0,10000,sx3event.Energy1,pcevent.Energy2*sinTheta_customV,pmlabel); + plotter->Fill2D("pmiscs_dPhi_sx3_PC"+tag,180,-360,360,180,-360,360,pcevent.pos.Phi()*180/M_PI,sx3event.pos.Phi()*180/M_PI,pmlabel); + plotter->Fill1D("pmiscs_dt_Anode_sx3_PC"+std::to_string(phicut)+tag,600,-2000,2000,pcevent.Time1-sx3event.Time1,pmlabel); + plotter->Fill1D("pmiscs_dt_Cathode_sx3"+tag,600,-2000,2000,pcevent.Time2-sx3event.Time1,pmlabel); + plotter->Fill2D("pmiscs_dt_Anode_E_sx3_PC"+std::to_string(phicut)+tag,600,-2000,2000,400,0,10,pcevent.Time1-sx3event.Time1,sx3event.Energy1,pmlabel); + plotter->Fill2D("pmiscs_dt_Cathode_E_sx3"+tag,600,-2000,2000,400,0,10,pcevent.Time2-sx3event.Time1,sx3event.Energy1,pmlabel); + plotter->Fill2D("pmiscs_dt_Cathodesx3_vsPCPhi"+tag,600,-2000,2000,180,-360,360,pcevent.Time2-sx3event.Time1,pcevent.pos.Phi()*180./M_PI,pmlabel); + plotter->Fill1D("pmiscs_pczfix"+tag,600,-300,300,pcz_fix,pmlabel); + plotter->Fill1D("pmiscs_pcz"+tag,600,-300,300,pcevent.pos.Z(),pmlabel); + + //double path_length_q = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; + //double path_length_s = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; + double path_length_s = (sx3event.pos-r_rhoMin_fix).Mag()*0.1; + double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1)-path_length_s); + + //plotter->Fill2D("sx3Ef_sx3E_matrix_all"+tag,400,0,10,400,0,10,sx3Efix,sx3event.Energy1,pmlabel); + plotter->Fill2D("pmiscs_dE3_Ef_Anodesx3"+tag,400,0,10,400,0,40000,sx3Efix,pcevent.Energy1*sinTheta_customV*3,pmlabel); + plotter->Fill2D("pmiscs_dE3_Ef_Cathodesx3"+tag,400,0,10,400,0,10000,sx3Efix,pcevent.Energy2*sinTheta_customV,pmlabel); + + plotter->Fill2D("pmiscs_Ef_vs_theta_sx3"+tag,100,0,180,800,0,20,theta_s*180/M_PI,sx3Efix,pmlabel); + plotter->Fill1D("pmiscs_VertexReconZ"+tag,800,-400,400,vertex_z,pmlabel); + plotter->Fill2D("pmiscs_VertexReconXY"+tag,200,-100,100,200,-100,100,r_rhoMin_fix.X(),r_rhoMin_fix.Y(),pmlabel); + plotter->Fill2D("pmiscs_VertexReconZ_vs_Ef"+tag,800,-400,400,800,0,20,vertex_z,sx3Efix,pmlabel); + plotter->Fill2D("pmiscs_VertexReconZ_vs_Ef"+tag+"_a"+std::to_string(pcevent.multi1),800,-400,400,800,0,20,vertex_z,sx3Efix,pmlabel); + }; + + plot_with_tag(); + if(cathode_alpha_select) + plot_with_tag("_cathode_alphas"); + else + plot_with_tag("_cathode_protons"); + + //plotter->Fill1D("pmisc_Ex_from_protons",200,-10,10,apkin_p.getExc(sx3Efix,theta_s*180/M_PI),pmlabel); + + }//end PCEvents loop + }//end sx3Events loop } @@ -1614,34 +1520,31 @@ void protonAlphaHistograms(HistPlotter* plotter, std::vector QQQ_Events, plotter->Fill1D("ap_qqq_sx3_dphi",180,-360,360,qqqevent.pos.Phi()*180/M_PI - sx3event.pos.Phi()*180/M_PI,aplabel); plotter->Fill2D("ap_qqq_sx3_dphi_vs_qqqphi",180,-360,360,180,-360,360,qqqevent.pos.Phi()*180/M_PI - sx3event.pos.Phi()*180/M_PI,qqqevent.pos.Phi()*180/M_PI,aplabel); plotter->Fill2D("ap_qqq_sx3_matrix",400,0,10,400,0,10,qqqevent.Energy1,sx3event.Energy1,aplabel); - + for(auto pcevent: PC_Events) { - - double pcz_fix = pcfix_func.Eval(pcevent.pos.Z())-15.0; + double pcz_fix = pcfix_func.Eval(pcevent.pos.Z())-5.0; TVector3 x2f(pcevent.pos.X(),pcevent.pos.Y(),pcz_fix); TVector3 x1(qqqevent.pos); 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; double vertex_z = r_rhoMin_fix.Z(); - double theta_q = (qqqevent.pos - TVector3(0,0,vertex_z)).Theta(); - //double theta_q = (qqqevent.pos - r_rhoMin_fix).Theta(); + //double theta_q = (qqqevent.pos - TVector3(0,0,vertex_z)).Theta(); + double theta_q = (qqqevent.pos - r_rhoMin_fix).Theta(); double sinTheta_customV = TMath::Sin(theta_q); - double theta_s = (sx3event.pos - TVector3(0,0,vertex_z)).Theta(); - //double theta_s = (sx3event.pos - r_rhoMin_fix).Theta(); + //double theta_s = (sx3event.pos - TVector3(0,0,vertex_z)).Theta(); + double theta_s = (sx3event.pos - r_rhoMin_fix).Theta(); double sinTheta_s = TMath::Sin(theta_s); //if(vertex_z<0 || vertex_z>100) continue; - - //double sinTheta = TMath::Sin((qqqevent.pos - pcevent.pos).Theta()); + + //double sinTheta = TMath::Sin((qqqevent.pos - pcevent.pos).Theta()); //plotter->Fill2D("sinTheta2_vs_sinTheta",80,-2,2,80,-2,2,sinTheta,sinTheta_customV,aplabel); - plotter->Fill2D("ap_dE_E_Anodesx3B",400,0,10,800,0,40000,sx3event.Energy1,pcevent.Energy1,aplabel); plotter->Fill2D("ap_dE_E_Cathodesx3B",400,0,10,800,0,10000,sx3event.Energy1,pcevent.Energy2,aplabel); plotter->Fill2D("ap_dE_E_AnodeQQQ",400,0,10,800,0,40000,qqqevent.Energy1,pcevent.Energy1,aplabel); plotter->Fill2D("ap_dE_E_CathodeQQQ",400,0,10,800,0,10000,qqqevent.Energy1,pcevent.Energy2,aplabel); plotter->Fill2D("ap_dE3_E_AnodeQQQ",400,0,10,400,0,40000,qqqevent.Energy1,pcevent.Energy1*sinTheta_customV,aplabel); - plotter->Fill2D("ap_dE3_E_CathodeQQQ",400,0,10,400,0,10000,qqqevent.Energy1,pcevent.Energy2*sinTheta_customV,aplabel); - + plotter->Fill2D("ap_dE3_E_CathodeQQQ",400,0,10,400,0,10000,qqqevent.Energy1,pcevent.Energy2*sinTheta_customV,aplabel); plotter->Fill2D("ap_dPhi_QQQ_PC",180,-360,360,180,-360,360,pcevent.pos.Phi()*180/M_PI,qqqevent.pos.Phi()*180/M_PI,aplabel); plotter->Fill2D("ap_dPhi_SX3_PC",180,-360,360,180,-360,360,pcevent.pos.Phi()*180/M_PI,sx3event.pos.Phi()*180/M_PI,aplabel); plotter->Fill1D("ap_dt_Anode_QQQ",600,-2000,2000,pcevent.Time1-qqqevent.Time1,aplabel); @@ -1650,20 +1553,20 @@ void protonAlphaHistograms(HistPlotter* plotter, std::vector QQQ_Events, plotter->Fill1D("ap_dt_Cathode_SX3",600,-2000,2000,pcevent.Time2-sx3event.Time1,aplabel); plotter->Fill1D("ap_pczfix",600,-300,300,pcz_fix,aplabel); plotter->Fill1D("ap_pcz",600,-300,300,pcevent.pos.Z(),aplabel); - - double path_length_q = (qqqevent.pos-TVector3(0,0,vertex_z)).Mag()*0.1; - double path_length_s = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; - //double path_length_q = (qqqevent.pos-r_rhoMin_fix).Mag()*0.1; - //double path_length_s = (sx3event.pos-r_rhoMin_fix).Mag()*0.1; - + + //double path_length_q = (qqqevent.pos-TVector3(0,0,vertex_z)).Mag()*0.1; + //double path_length_s = (sx3event.pos-TVector3(0,0,vertex_z)).Mag()*0.1; + double path_length_q = (qqqevent.pos-r_rhoMin_fix).Mag()*0.1; + double path_length_s = (sx3event.pos-r_rhoMin_fix).Mag()*0.1; //We know that alphas predominantly are detected in QQQs, and protons in SX3s, and that protons don't leave much of a trace in dE layer. //Using the estimated path lengths, we correct alpha eloss in qqq, and protons in sx3. The result should (hopefully be) vertex independent. - double qqqEfix = cm_to_MeV->Eval(MeV_to_cm->Eval(qqqevent.Energy1)-path_length_q); double sx3Efix = cm_to_MeVp->Eval(MeV_to_cm_p->Eval(sx3event.Energy1)-path_length_s); //plotter->Fill2D("qqqEf_sx3E_matrix_all",400,0,10,400,0,10,qqqEfix,sx3event.Energy1,aplabel); + plotter->Fill2D("ap_dE3_Ef_AnodeQQQ",400,0,10,400,0,40000,qqqEfix,pcevent.Energy1*sinTheta_customV,aplabel); + plotter->Fill2D("ap_dE3_Ef_CathodeQQQ",400,0,10,400,0,10000,qqqEfix,pcevent.Energy2*sinTheta_customV,aplabel); + plotter->Fill2D("ap_qqqEf_sx3Ef_matrix",400,0,10,400,0,10,qqqEfix,sx3Efix,aplabel); - plotter->Fill2D("ap_Ef_vs_theta_qqq",100,0,180,400,0,10,theta_q*180/M_PI,qqqEfix,aplabel); plotter->Fill2D("ap_Ef_vs_theta_sx3",100,0,180,400,0,10,theta_s*180/M_PI,sx3Efix,aplabel); plotter->Fill2D("ap_theta_vs_theta_qqq_sx3",100,0,180,100,0,180,theta_q*180/M_PI,theta_s*180/M_PI,aplabel); @@ -1671,7 +1574,6 @@ void protonAlphaHistograms(HistPlotter* plotter, std::vector QQQ_Events, plotter->Fill2D("ap_VertexReconXY",200,-100,100,200,-100,100,r_rhoMin_fix.X(),r_rhoMin_fix.Y(),aplabel); plotter->Fill1D("ap_Ex_from_protons",200,-10,10,apkin_p.getExc(sx3Efix,theta_s*180/M_PI),aplabel); plotter->Fill1D("ap_Ex_from_alpha",200,-10,10,apkin_a.getExc(qqqEfix,theta_q*180/M_PI),aplabel); - if(pcevent.multi1==1 && pcevent.multi2==2) { //one-anode, two-cathode events, as originally intended //std::cout << "Test" << std::endl; plotter->Fill1D("ap_VertexReconZ_a1c2",400,-200,200,vertex_z,aplabel); @@ -1680,24 +1582,19 @@ void protonAlphaHistograms(HistPlotter* plotter, std::vector QQQ_Events, plotter->Fill2D("ap_Ef_vs_theta_qqq_a1c2",100,0,180,400,0,10,theta_q*180/M_PI,qqqEfix,aplabel); plotter->Fill1D("ap_Ex_from_protons_a1c2",200,-10,10,apkin_p.getExc(sx3Efix,theta_s*180/M_PI),aplabel); plotter->Fill1D("ap_Ex_from_alpha_a1c2",200,-10,10,apkin_a.getExc(qqqEfix,theta_q*180/M_PI),aplabel); - - //std::cout << apkin_p.getExc(sx3Efix,theta_s*180/M_PI) << " " << apkin_a.getExc(qqqEfix,theta_q*180/M_PI)<< std::endl; + + //std::cout << apkin_p.getExc(sx3Efix,theta_s*180/M_PI) << " " << apkin_a.getExc(qqqEfix,theta_q*180/M_PI)<< std::endl; plotter->Fill2D("ap_Ef_vs_theta_sx3_a1c2",100,0,180,400,0,10,theta_s*180/M_PI,sx3Efix,aplabel); - //plotter->Fill2D("qqqEf_sx3E_matrix",400,0,10,400,0,10,qqqEfix,sx3event.Energy1,aplabel); plotter->Fill2D("ap_qqq_sx3_matrix_a1c2",400,0,10,400,0,10,qqqevent.Energy1,sx3event.Energy1,aplabel); plotter->Fill2D("ap_qqqEf_sx3Ef_matrix_a1c2",400,0,10,400,0,10,qqqEfix,sx3Efix,aplabel); //std::cout << sx3event.Energy1 << " " << path_length_s << " " << sx3Efix << std::endl; - + //plotter->Fill2D("dE3_Ef_AnodeQQQ_a1c2",400,0,10,400,0,40000,qqqEfix,pcevent.Energy1*sinTheta_customV,aplabel); //plotter->Fill2D("dE3_Ef_CathodeQQQ_a1c2",400,0,10,400,0,10000,qqqEfix,pcevent.Energy2*sinTheta_customV,aplabel); - } //end if(a1c2) loop }//end PC_Events for loop - - }//end SX3_Events for loop } //end QQQ_Events for loop, end sidetrack a(p,p) return; } - diff --git a/MakeVertex.h b/MakeVertex.h index d14bf2e..5e4adae 100644 --- a/MakeVertex.h +++ b/MakeVertex.h @@ -124,11 +124,6 @@ void MakeVertex::Init(TTree *tree) { fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh); fChain->SetBranchAddress("pcE", &pc.e, &b_pcE); fChain->SetBranchAddress("pcT", &pc.t, &b_pcT); - /*fChain->SetBranchAddress("miscMulti", &misc.multi, &b_miscMulti); - fChain->SetBranchAddress("miscID", &misc.id, &b_miscID); - fChain->SetBranchAddress("miscCh", &misc.ch, &b_miscCh); - fChain->SetBranchAddress("miscE", &misc.e, &b_miscE); - fChain->SetBranchAddress("miscT", &misc.t, &b_miscT);*/ } Bool_t MakeVertex::Notify() { diff --git a/cutlist.txt b/cutlist.txt new file mode 100644 index 0000000..ee225c9 --- /dev/null +++ b/cutlist.txt @@ -0,0 +1,4 @@ +results/anode_strip_diagonalCut1.root anode_sx3_diag1 +results/anode_strip_diagonalCut2.root anode_sx3_diag2 +results/anode_wedge_diagonalCut1.root anode_qqq_diag1 +results/anode_wedge_diagonalCut2.root anode_qqq_diag2 diff --git a/results/compare.C b/results/compare.C index 04ddf14..597a3e8 100644 --- a/results/compare.C +++ b/results/compare.C @@ -5,7 +5,7 @@ //TFile fin("4He(pp)_candidate_kinematic_curve.root"); //TH2F *h2 = (TH2F*)(fin.Get("Ef_thetaf_AnodeQQQR_TC1_PC1_pidlow0_1")); -// TFile fin("../results_run15.root"); + //TFile fin("../results_run017.root"); TFile fin("out.root"); //TH2F *h2 = (TH2F*)(fin.Get("Ef_thetaf_AnodeQQQR_TC1_PC1_pidlow0_1")); // TFile fin("run15_sx3_qqq_alphas_kinematic.root"); @@ -20,7 +20,7 @@ TGraph angles("a(p,p)a_kinematics_7MeV_p.txt","%lf %*lf %lf %*lf"); TGraph angles2(angles.GetN(),angles.GetY(),angles.GetX()); h2 = (TH2F*)(fin.Get("a(p,p)/ap_theta_vs_theta_qqq_sx3_a1c2")); - h2->Draw(); + h2->Draw("colz"); angles2.Draw("L SAME"); gPad->Modified(); gPad->Update(); diff --git a/results/compareE.C b/results/compareE.C index f3f7928..b55f126 100644 --- a/results/compareE.C +++ b/results/compareE.C @@ -6,7 +6,6 @@ //TFile fin("4He(pp)_candidate_kinematic_curve.root"); //TH2F *h2 = (TH2F*)(fin.Get("Ef_thetaf_AnodeQQQR_TC1_PC1_pidlow0_1")); -// TFile fin("../results_run18.root"); // TFile fin("out_zminus5_catima.root"); TFile fin("out.root"); //TH2F *h2 = (TH2F*)(fin.Get("Ef_thetaf_AnodeQQQR_TC1_PC1_pidlow0_1")); diff --git a/results/compareE_pcalpha.C b/results/compareE_pcalpha.C new file mode 100644 index 0000000..e7a7cc2 --- /dev/null +++ b/results/compareE_pcalpha.C @@ -0,0 +1,26 @@ +{ + TCanvas c; + //c.SetLogy(kTRUE); + TGraph kinematics("a(p,p)a_kinematics_7MeV_p.txt","%*lf %*lf %lf %lf"); + kinematics.Scale(4.0); + kinematics.SetLineColor(kBlue); + kinematics.SetLineWidth(2.); + //TFile fin("4He(pp)_candidate_kinematic_curve.root"); + //TH2F *h2 = (TH2F*)(fin.Get("Ef_thetaf_AnodeQQQR_TC1_PC1_pidlow0_1")); + +// TFile fin("../results_run017.root"); +// TFile fin("out_zminus5_catima.root"); + TFile fin("out.root"); + //TH2F *h2 = (TH2F*)(fin.Get("Ef_thetaf_AnodeQQQR_TC1_PC1_pidlow0_1")); +// TFile fin("run15_sx3_qqq_alphas_kinematic.root"); + TH2F *h2 = (TH2F*)(fin.Get("proton+misc_cathode_alphas/pmisc_Ef_vs_theta_qqq_cathode_alphas")); + h2->Draw("col"); + h2->RebinY(8); + kinematics.Draw("L SAME"); + gPad->Modified(); + gPad->Update(); + while(gPad->WaitPrimitive()); + + gPad->SaveAs("c1.png"); + +} diff --git a/results/compareE_pcalpha_ow.C b/results/compareE_pcalpha_ow.C new file mode 100644 index 0000000..e140f24 --- /dev/null +++ b/results/compareE_pcalpha_ow.C @@ -0,0 +1,26 @@ +{ + TCanvas c; + //c.SetLogy(kTRUE); + TGraph kinematics("a(p,p)a_kinematics_7MeV_p.txt","%*lf %*lf %lf %lf"); + kinematics.Scale(4.0); + kinematics.SetLineColor(kBlue); + kinematics.SetLineWidth(2.); + //TFile fin("4He(pp)_candidate_kinematic_curve.root"); + //TH2F *h2 = (TH2F*)(fin.Get("Ef_thetaf_AnodeQQQR_TC1_PC1_pidlow0_1")); + +// TFile fin("../results_run017.root"); +// TFile fin("out_zminus5_catima.root"); + TFile fin("out.root"); + //TH2F *h2 = (TH2F*)(fin.Get("Ef_thetaf_AnodeQQQR_TC1_PC1_pidlow0_1")); +// TFile fin("run15_sx3_qqq_alphas_kinematic.root"); + TH2F *h2 = (TH2F*)(fin.Get("ainterp_noc/pmisc_ow_Ef_vs_theta_qqq")); + h2->Draw("col"); + h2->RebinY(8); + kinematics.Draw("L SAME"); + gPad->Modified(); + gPad->Update(); + while(gPad->WaitPrimitive()); + + gPad->SaveAs("c1.png"); + +} diff --git a/results/haddnow b/results/haddnow new file mode 100644 index 0000000..49fc057 --- /dev/null +++ b/results/haddnow @@ -0,0 +1,6 @@ +rm out.root +#hadd out.root ../results_run15.root ../results_run18.root ../results_run19.root ../results_run17.root ../results_run20.root ../results_run21.root ../results_run22.root +#hadd out.root ../results_run015.root ../results_run018.root ../results_run019.root ../results_run017.root ../results_run020.root ../results_run021.root ../results_run022.root +hadd out.root ../results_run018.root ../results_run019.root ../results_run017.root ../results_run020.root ../results_run021.root ../results_run022.root + +#hadd out.root ../results_run18.root ../results_run19.root ../results_run17.root ../results_run20.root ../results_run21.root ../results_run22.root diff --git a/run_17F.sh b/run_17F.sh index 906fb88..7fec7ef 100644 --- a/run_17F.sh +++ b/run_17F.sh @@ -1,46 +1,34 @@ rm results_run*.root export DATASET="17F" -export flip180="0" -export flipa=0 -export anode_offset=0 -#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; -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_007_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run07.root; -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_008_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run08.root; -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_009_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run09.root; -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_010_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run10.root; -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_011_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run11.root; -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_012_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run12.root; -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_013_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run13.root; -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_014_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run14.root; - -#17F pulser runs -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/PulserRun_015_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run15.root; -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/PulserRun_016_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run16.root; -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/PulserRun_017_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run17.root; - -#17F alpha run with gas -#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_018_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run18.root; -#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; -#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; -#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; +export reactiondata=1 #17F reaction data -export flip180="0" -declare -i run=231 #49 -while [[ $run -lt 258 ]]; do #392 - wrun=$(printf "%03d" $run) - file_exists=$(test -f ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root) - if [[ $file_exists -ne 0 ]]; then continue; fi - root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run$wrun.root; - run=run+1 -done +#declare -i run=231 #49 +#while [[ $run -lt 258 ]]; do #392 +# wrun=$(printf "%03d" $run) +# file_exists=$(test -f ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root) +# if [[ $file_exists -ne 0 ]]; then continue; fi +# root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run$wrun.root; +# run=run+1 +#done +function run_once() { + wrun=$(printf "%03d" $1) + file_exists=$(test -f /home/sud/Desktop/Software2/ANASEN_analysis/data/17F_fsu_files/Run_"$wrun"_mapped.root) + if [[ $file_exists -ne 0 ]]; then return; fi + root -q -l -b -x /home/sud/Desktop/Software2/ANASEN_analysis/data/17F_fsu_files/Run_"$wrun"_mapped.root -e $(printf 'tree->Process("MakeVertex.C+O","analyzed_run%s.root")' "$wrun"); + mv analyzed_run$wrun.root results_run$wrun.root; +} + +export -f run_once +#run_once 351 +parallel -j 6 --ctag run_once {1} ::: {350..400} rm output_17F.root -hadd -j 4 -k output_17F.root results_run2*.root +hadd -j 4 -k output_17F.root results_run3*.root unset souce_vertex unset DATASET unset flip180 unset flipa unset anode_offset +unset reactiondata diff --git a/run_27Al.sh b/run_27Al.sh index 47b8937..b669a96 100644 --- a/run_27Al.sh +++ b/run_27Al.sh @@ -2,7 +2,10 @@ export DATASET="27Al" export flip180="0" export flipa=0 +export reactiondata=1 export anode_offset=0 +root -l -q -x -e ".L MakeVertex.C++" +rm 27Al_output/*.root #declare -i run=28 #while [[ $run -lt 34 ]]; do #runs 1 to 84 # wrun=$(printf "%03d" $run) @@ -10,21 +13,29 @@ export anode_offset=0 # run=run+1 #done -declare -i run=50 -while [[ $run -lt 59 ]]; do #runs 1 to 84 +declare -i run=78 +while [[ $run -lt 89 ]]; do #runs 1 to 84 wrun=$(printf "%03d" $run) - root -q -l -b -x ../ANASEN_analysis/data/27Al_Data/Run_"$wrun"_mapped.root -e 'tree->Process("MakeVertex.C+O","Analyzer_27Al.root")'; mv Analyzer_27Al.root 27Al_output/results_run$wrun.root; + #root -q -l -b -x ../ANASEN_analysis/data/27Al_Data/Run_"$wrun"_mapped.root -e 'tree->Process("MakeVertex.C+O","Analyzer_27Al.root")'; mv Analyzer_27Al.root 27Al_output/results_run$wrun.root; run=run+1 done -rm output.root +function run_once() { + wrun=$(printf "%02d" $1) + file_exists=$(test -f ../ANASEN_analysis/data/27Al_Data/Run_0"$wrun"_mapped.root) + if [[ $file_exists -ne 0 ]]; then return; fi + root -q -l -b -x ../ANASEN_analysis/data/27Al_Data/Run_0"$wrun"_mapped.root -e $(printf 'tree->Process("MakeVertex.C+O","analyzed_run%s.root")' "$wrun"); + mv analyzed_run$wrun.root 27Al_output/results_run$wrun.root; +} + +export -f run_once +time parallel -j 6 --ctag run_once {1} ::: {50..59} +time parallel -j 2 --ctag run_once {1} ::: {78..89} + + hadd -k -j 4 output.root 27Al_output/results_run*.root mv output.root output_27Al.root -#root -q -l -b -x -e '.L MakeVertex.C+O'; -#halfproc=3 -#parallel --ctag --bar -j $halfproc ./run.sh ::: {028..034} ::: 27Al #color-tag, linebuffer, then run run.sh in parallel - unset souce_vertex unset DATASET unset flip180 diff --git a/run_sx3.sh b/run_sx3.sh index 60d64ab..914edf5 100644 --- a/run_sx3.sh +++ b/run_sx3.sh @@ -1,8 +1,9 @@ #Alpha runs at different spacer positions #rm results_run*.root -export flipa=0 export anode_offset=0 +export cathode_offset=0 export DATASET="27Al" +root -l -q -x -e ".L MakeVertex.C++" if [[ 1 -eq 0 ]]; then #root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run09.root; root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_001_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run01.root; @@ -18,9 +19,8 @@ fi #exit #alpha+gas 27Al 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 1 ]]; then +if [[ 1 -eq 0 ]]; then #export timecut_low=230.0; export timecut_low=400.0; #export timecut_high=400.0; @@ -31,18 +31,15 @@ export timecut_low=400.0; 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 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; +#export source_vertex=-57.28; 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; +export source_vertex=-135.68; 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; +exit 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; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_020_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run20.root; @@ -51,6 +48,18 @@ root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_022_mapped.root -e 'tree->Pr exit fi +if [[ 1 -eq 1 ]]; then +function run_once() { + wrun=$(printf "%03d" $1) + file_exists=$(test -f ../ANASEN_analysis/data/27Al_Data/Run_"$wrun"_mapped.root) + if [[ $file_exists -ne 0 ]]; then return; fi + root -q -l -b -x ../ANASEN_analysis/data/27Al_Data/Run_"$wrun"_mapped.root -e $(printf 'tree->Process("MakeVertex.C+O","analyzed_run%s.root")' "$wrun"); + mv analyzed_run$wrun.root results_run$wrun.root; +} +export -f run_once +time parallel -j 4 --ctag run_once {1} ::: {16..22} +exit +fi #27Al reaction data #root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_051_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run51.root; #root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_078_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run78.root; @@ -59,7 +68,6 @@ fi #root -l -x results_run19.root results_run12.root -e "new TBrowser" #exit export DATASET="17F" -#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; @@ -86,7 +94,6 @@ export source_vertex=-73.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/S exit fi #17F reaction data -#export flip180="0" if [[ 1 -eq 0 ]]; then export source_vertex=-57.28; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/ProtonRun_035_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run35.root; #export source_vertex=-8.28; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/ProtonRun_036_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root resulrs_run36.root; @@ -101,11 +108,8 @@ export source_vertex=-57.28; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/P #root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Run_104_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run104.root; #mv Analyzer_SX3.root results_run19.root; fi -unset flipa -unset flipc unset anode_offset unset cathode_offset unset souce_vertex unset DATASET -unset flip180 unset timecut_low, timecut_high diff --git a/scratch/sx3z_vs_phiz/compare_wire_offsets.C b/scratch/sx3z_vs_phiz/compare_wire_offsets.C new file mode 100644 index 0000000..f397bc1 --- /dev/null +++ b/scratch/sx3z_vs_phiz/compare_wire_offsets.C @@ -0,0 +1,107 @@ +#include "testmodel.h" + +int quit=0; +void handler(int){quit=0;} + +int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5}; +void compare_wire_offsets(){ + signal(SIGINT,handler); + TCanvas c("c1","c1",0,0,1600,800); + c.Divide(2,1); + + TF1 f1("model",model,-200,200,2); + f1.SetNpx(10000); + std::vector pars = {0.0,1.}; + f1.SetParameters(pars.data()); + f1.SetLineColor(kGreen+2); + f1.SetLineStyle(kLine); + + + + + TFile* f=NULL; + std::vector files; + int ctr=0; + for(int i=12; i<=21; i++) { + 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")); + std::cout << "aaa" << h23 << std::endl; + h23->SetLineColorAlpha(kOrange,0.75); + 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)); + //h2->Draw("colz same"); + h2->SetLineColorAlpha(colors[ctr],0.75); + h2->Draw("box same"); + } + TF1 eqline("x","x",-200,200); + //eqline.Draw("SAME"); + //f1.Draw("same"); + + c1->Modified(); + c1->Update(); + + + auto c2=c.cd(2); + c2->SetGrid(1,1); + + /*TH2F *h3 = (TH2F*)(f->Get("sx3phi_vs_pcphi1")); + if(!h3) continue; + h3->SetTitle(Form("case%d",i)); + h3->Draw("colz"); + eqline.Draw("SAME");*/ + + if(i==12) { + //TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int")); + TH2F *h23 = (TH2F*)(f->Get("pczfix_vs_qqqpczguess_A1C2")); + std::cout << "aaa" << h23 << std::endl; + h23->SetLineColorAlpha(kOrange,0.75); + 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("pczfix_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)); + //h2->Draw("colz same"); + h2->SetLineColorAlpha(colors[ctr],0.75); + h2->Draw("box same"); + } + + + c2->Modified(); + c2->Update(); + ctr+=1; + + while(gPad->WaitPrimitive()); + + files.emplace_back(f); + std::cout <<"Test" << std::endl; + if(i==21) { + i=11; + c.Modified(); c.Update(); + c.SaveAs("test.png"); + c.Clear(); + c.Divide(2,1); + ctr=0; + } + //if(quit) break; + } + for(auto file : files) { + file->Close(); + } +} diff --git a/scratch/sx3z_vs_phiz/compare_wire_offsets_c3.C b/scratch/sx3z_vs_phiz/compare_wire_offsets_c3.C new file mode 100644 index 0000000..5af0468 --- /dev/null +++ b/scratch/sx3z_vs_phiz/compare_wire_offsets_c3.C @@ -0,0 +1,88 @@ +#include "testmodel.h" + +int quit=0; +void handler(int){quit=0;} + +int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5}; +void compare_wire_offsets_c3(){ + signal(SIGINT,handler); + TCanvas c("c1","c1",0,0,1600,800); + c.Divide(2,1); + + TF1 f1("model",model,-200,200,2); + f1.SetNpx(10000); + std::vector pars = {0.0,1.}; + f1.SetParameters(pars.data()); + f1.SetLineColor(kGreen+2); + f1.SetLineStyle(kLine); + + + + + TFile* f=NULL; + std::vector files; + int ctr=0; + for(int i=12; i<=21; i++) { + 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_A1C3")); + if(h23) { + std::cout << "aaa" << h23 << std::endl; + h23->SetLineColorAlpha(kOrange,0.75); + 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_A1C3")); + //TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(h2) { + h2->SetTitle(Form("case%d",i)); + //h2->Draw("colz same"); + h2->SetLineColorAlpha(colors[ctr],0.75); + h2->Draw("box same"); + } + } + TF1 eqline("x","x",-200,200); + //eqline.Draw("SAME"); + //f1.Draw("same"); + + c1->Modified(); + c1->Update(); + + + auto c2=c.cd(2); + c2->SetGrid(1,1); + + /*TH2F *h3 = (TH2F*)(f->Get("sx3phi_vs_pcphi1")); + if(!h3) continue; + h3->SetTitle(Form("case%d",i)); + h3->Draw("colz"); + eqline.Draw("SAME");*/ + + c2->Modified(); + c2->Update(); + ctr+=1; + + while(gPad->WaitPrimitive()); + + files.emplace_back(f); + std::cout <<"Test" << std::endl; + if(i==21) { + i=11; + c.Modified(); c.Update(); + c.SaveAs("test.png"); + c.Clear(); + c.Divide(2,1); + ctr=0; + } + //if(quit) break; + } + for(auto file : files) { + file->Close(); + } +} diff --git a/scratch/sx3z_vs_phiz/scan_offset.C b/scratch/sx3z_vs_phiz/scan_offset.C index f6ea1b0..eb7ba91 100644 --- a/scratch/sx3z_vs_phiz/scan_offset.C +++ b/scratch/sx3z_vs_phiz/scan_offset.C @@ -1,72 +1,93 @@ -#include "testmodel.h" - int quit=0; void handler(int){quit=0;} - int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5}; void scan_offset(){ signal(SIGINT,handler); TCanvas c("c1","c1",0,0,1600,800); - c.Divide(2,1); - - TF1 f1("model",model,-200,200,2); - f1.SetNpx(10000); - std::vector pars = {0.0,1.}; - f1.SetParameters(pars.data()); - f1.SetLineColor(kGreen+2); - f1.SetLineStyle(kLine); - - - + c.Divide(2,2); TFile* f=NULL; std::vector files; int ctr=0; for(int i=12; i<=21; i++) { + if(i>=13 && i<=17) continue; 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")); + TH2F *h23 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int")); + //TH2F *h23 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int_self")); std::cout << "aaa" << h23 << std::endl; h23->SetLineColorAlpha(kOrange,0.75); h23->GetYaxis()->SetRangeUser(-200,200); - h23->Draw("box"); - while(gPad->WaitPrimitive()); - } else { + h23->Draw("col"); + } { //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")); + //TH2F *h2 = (TH2F*)(f->Get("pcz_vs_sx3pczguess_self")); std::cout << h2 << std::endl; //TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); - if(!h2) continue; - h2->SetTitle(Form("case%d",i)); - //h2->Draw("colz same"); - h2->SetLineColorAlpha(colors[ctr],0.75); - h2->Draw("box same"); - f1.Draw("same"); + if(h2) { + h2->SetTitle(Form("case%d",i)); + //h2->Draw("colz same"); + h2->SetLineColorAlpha(colors[ctr],0.75); + h2->Draw("col same"); + } } - TF1 eqline("x","x",-200,200); - eqline.Draw("SAME"); - c1->Modified(); c1->Update(); ctr+=1; + auto c1a = c.cd(2); + c1a->SetGrid(1,1); + TH2F *h24 = (TH2F*)(f->Get("sx3phi_vs_pcphi1")); + if(h24) { + h24->Draw("box"); + h24->SetLineColor(kRed); + } + TH2F *h25 = (TH2F*)(f->Get("phiPC_vs_phiQQQ_TimeCut")); + if(h25) { + h25->SetLineColor(kBlue); + h25->Draw("box same"); + } + c1a->Modified(); + c1a->Update(); - auto c2=c.cd(2); + + auto c2=c.cd(3); c2->SetGrid(1,1); - TH2F *h3 = (TH2F*)(f->Get("sx3phi_vs_pcphi1")); + TH2F *h3 = (TH2F*)(f->Get("d_sx3pczguess_minus_pcz_a1c2")); // TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); - if(!h3) continue; - h3->SetTitle(Form("case%d",i)); - h3->Draw("colz"); - eqline.Draw("SAME"); + if(h3) { + h3->SetLineColor(kRed); + h3->SetTitle(Form("z_vs_guess_run%d",i)); + h3->Draw("colz"); + } + TH2F *h4 = (TH2F*)(f->Get("d_qqqpczguess_minus_pcz_a1c2")); + //TH2F *h4 = (TH2F*)(f->Get("d_sx3pczguess_minus_pczfix_a1c2")); + if(h4) + h4->Draw("colz same"); c2->Modified(); c2->Update(); + auto c23 = c.cd(4); + c23->SetGrid(1,1); + TH1F *hdiff2 = (TH1F*)(f->Get("d_sx3phi_minus_pcphi1")); + if(hdiff2) { + hdiff2->SetLineColor(kRed); + hdiff2->Draw(""); + } + TH1F *hdiff = (TH1F*)(f->Get("d_phiPC_phiQQQ_TimeCut")); + if(hdiff) { + hdiff->SetLineColor(kBlue); + hdiff->Draw("SAME"); + } + + c23->Modified(); + c23->Update(); while(gPad->WaitPrimitive()); files.emplace_back(f); @@ -74,7 +95,7 @@ void scan_offset(){ if(i==21) { i=11; c.Clear(); - c.Divide(2,1); + c.Divide(2,2); ctr=0; } //if(quit) break; diff --git a/scratch/sx3z_vs_phiz/scan_offset_a2c1.C b/scratch/sx3z_vs_phiz/scan_offset_a2c1.C new file mode 100644 index 0000000..2a1807d --- /dev/null +++ b/scratch/sx3z_vs_phiz/scan_offset_a2c1.C @@ -0,0 +1,86 @@ +#include "testmodel.h" + +int quit=0; +void handler(int){quit=0;} + +int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5}; +void scan_offset_a2c1(){ + signal(SIGINT,handler); + TCanvas c("c1","c1",0,0,1600,800); + c.Divide(2,1); + + TF1 f1("model",model,-200,200,2); + f1.SetNpx(10000); + std::vector pars = {0.0,1.}; + f1.SetParameters(pars.data()); + f1.SetLineColor(kGreen+2); + f1.SetLineStyle(kLine); + + + + + TFile* f=NULL; + std::vector files; + int ctr=0; + for(int i=12; i<=21; i++) { + 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("pcz_vs_qqqpczguess_A2C1")); + std::cout << "aaa" << h23 << std::endl; + h23->SetLineColorAlpha(kOrange,0.75); + h23->GetYaxis()->SetRangeUser(-200,200); + h23->Draw("box"); + while(gPad->WaitPrimitive()); + } + { + //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_A2C1")); + std::cout << h2 << std::endl; + //TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(!h2) continue; + h2->SetTitle(Form("case%d",i)); + //h2->Draw("colz same"); + h2->SetLineColorAlpha(colors[ctr],0.75); + h2->Draw("box same"); + f1.Draw("same"); + } + TF1 eqline("x","x",-200,200); + eqline.Draw("SAME"); + + c1->Modified(); + c1->Update(); + ctr+=1; + + + auto c2=c.cd(2); + c2->SetGrid(1,1); + + TH2F *h3 = (TH2F*)(f->Get("sx3phi_vs_pcphi1")); +// TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(!h3) continue; + h3->SetTitle(Form("case%d",i)); + h3->Draw("colz"); + eqline.Draw("SAME"); + c2->Modified(); + c2->Update(); + + 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; + } + for(auto file : files) { + file->Close(); + } +} diff --git a/scratch/sx3z_vs_phiz/scan_offset_ainterp.C b/scratch/sx3z_vs_phiz/scan_offset_ainterp.C new file mode 100644 index 0000000..31de54c --- /dev/null +++ b/scratch/sx3z_vs_phiz/scan_offset_ainterp.C @@ -0,0 +1,109 @@ +int quit=0; +void handler(int){quit=0;} +int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5}; +void scan_offset_ainterp(){ + signal(SIGINT,handler); + TCanvas c("c1","c1",0,0,1600,800); + c.Divide(2,2); + + TFile* f=NULL; + std::vector files; + int ctr=0; + for(int i=12; i<=21; i++) { + if(i>=13 && i<=17) continue; + 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("ainterp_noc/pcZ_ainterp_sx3pczguess_TC1_ignC2")); + //TH2F *h23 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int_self")); + std::cout << "aaa" << h23 << std::endl; + if(h23) { + h23->SetLineColorAlpha(kOrange,0.75); + h23->GetYaxis()->SetRangeUser(-200,200); + h23->Draw("col"); + } + } + { + //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("ainterp_noc/pcZ_ainterp_sx3pczguess_TC1_ignC2")); + //TH2F *h2 = (TH2F*)(f->Get("pcz_vs_sx3pczguess_self")); + std::cout << h2 << std::endl; + //TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(h2) { + h2->SetTitle(Form("case%d",i)); + //h2->Draw("colz same"); + h2->SetLineColorAlpha(colors[ctr],0.75); + h2->Draw("col same"); + } + } + c1->Modified(); + c1->Update(); + ctr+=1; + + auto c1a = c.cd(2); + c1a->SetGrid(1,1); + TH2F *h24 = (TH2F*)(f->Get("sx3phi_vs_pcphi1")); + if(h24) { + h24->Draw("box"); + h24->SetLineColor(kRed); + } + TH2F *h25 = (TH2F*)(f->Get("phiPC_vs_phiQQQ_TimeCut")); + if(h25) { + h25->SetLineColor(kBlue); + h25->Draw("box same"); + } + c1a->Modified(); + c1a->Update(); + + + auto c2=c.cd(3); + c2->SetGrid(1,1); + + TH2F *h3 = (TH2F*)(f->Get("d_sx3pczguess_minus_pcz_a1c2")); +// TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(h3) { + h3->SetLineColor(kRed); + h3->SetTitle(Form("z_vs_guess_run%d",i)); + h3->Draw("colz"); + } + TH2F *h4 = (TH2F*)(f->Get("d_qqqpczguess_minus_pcz_a1c2")); + //TH2F *h4 = (TH2F*)(f->Get("d_sx3pczguess_minus_pczfix_a1c2")); + if(h4) + h4->Draw("colz same"); + c2->Modified(); + c2->Update(); + + auto c23 = c.cd(4); + c23->SetGrid(1,1); + TH1F *hdiff2 = (TH1F*)(f->Get("d_sx3phi_minus_pcphi1")); + if(hdiff2) { + hdiff2->SetLineColor(kRed); + hdiff2->Draw(""); + } + TH1F *hdiff = (TH1F*)(f->Get("d_phiPC_phiQQQ_TimeCut")); + if(hdiff) { + hdiff->SetLineColor(kBlue); + hdiff->Draw("SAME"); + } + + c23->Modified(); + c23->Update(); + while(gPad->WaitPrimitive()); + + files.emplace_back(f); + std::cout <<"Test" << std::endl; + if(i==21) { + i=11; + c.Clear(); + c.Divide(2,2); + ctr=0; + } + //if(quit) break; + } + for(auto file : files) { + file->Close(); + } +} diff --git a/scratch/sx3z_vs_phiz/scan_offset_fix.C b/scratch/sx3z_vs_phiz/scan_offset_fix.C index 20dde89..a1803e6 100644 --- a/scratch/sx3z_vs_phiz/scan_offset_fix.C +++ b/scratch/sx3z_vs_phiz/scan_offset_fix.C @@ -30,10 +30,13 @@ void scan_offset_fix(){ if(i==12) { //TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int")); TH2F *h23 = (TH2F*)(f->Get("pczfix_vs_qqqpczguess_A1C2")); + h23->RebinX(2); h23->SetLineColorAlpha(kOrange,0.75); - h23->Draw("box SAME"); + h23->GetYaxis()->SetRangeUser(-80,120); + h23->GetXaxis()->SetRangeUser(-80,120); + h23->Draw("col SAME"); - } 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("pczfix_vs_sx3pczguess_A1C2")); @@ -42,7 +45,7 @@ void scan_offset_fix(){ h2->SetTitle(Form("case%d",i)); //h2->Draw("colz same"); h2->SetLineColorAlpha(colors[ctr],0.75); - h2->Draw("box same"); + h2->Draw("col same"); //f1.Draw("same"); } TF1 eqline("x","x",-200,200); diff --git a/scratch/sx3z_vs_phiz/scan_offset_fix_1d.C b/scratch/sx3z_vs_phiz/scan_offset_fix_1d.C new file mode 100644 index 0000000..a150252 --- /dev/null +++ b/scratch/sx3z_vs_phiz/scan_offset_fix_1d.C @@ -0,0 +1,105 @@ +int quit=0; +void handler(int){quit=0;} +int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5}; +void scan_offset_fix_1d(){ + signal(SIGINT,handler); + TCanvas c("c1","c1",0,0,1600,800); + c.Divide(2,2); + + TFile* f=NULL; + std::vector files; + int ctr=0; + for(int i=12; i<=21; i++) { + if(i>=13 && i<=17) continue; + 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("pczfix_vs_qqqpczguess_A1C2")); + //TH2F *h23 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int_self")); + std::cout << "aaa" << h23 << std::endl; + h23->SetLineColorAlpha(kOrange,0.75); + h23->GetYaxis()->SetRangeUser(-200,200); + h23->Draw("col"); + } { + //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("pczfix_vs_sx3pczguess_A1C2")); + //TH2F *h2 = (TH2F*)(f->Get("pcz_vs_sx3pczguess_self")); + std::cout << h2 << std::endl; + //TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(h2) { + h2->SetTitle(Form("case%d",i)); + //h2->Draw("colz same"); + h2->SetLineColorAlpha(colors[ctr],0.75); + h2->Draw("col same"); + } + } + c1->Modified(); + c1->Update(); + ctr+=1; + + auto c1a = c.cd(2); + c1a->SetGrid(1,1); + TH2F *h24 = (TH2F*)(f->Get("sx3phi_vs_pcphi1")); + if(h24) { + h24->Draw("box"); + h24->SetLineColor(kRed); + } + TH2F *h25 = (TH2F*)(f->Get("phiPC_vs_phiQQQ_TimeCut")); + if(h25) { + h25->SetLineColor(kBlue); + h25->Draw("box same"); + } + c1a->Modified(); + c1a->Update(); + + + auto c2=c.cd(3); + c2->SetGrid(1,1); + + TH2F *h3 = (TH2F*)(f->Get("d_sx3pczguess_minus_pczfix_a1c2")); +// TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(h3) { + h3->SetLineColor(kRed); + h3->SetTitle(Form("z_vs_guess_run%d",i)); + h3->Draw("colz"); + } + TH2F *h4 = (TH2F*)(f->Get("d_qqqpczguess_minus_pczfix_a1c2")); + if(h4) + h4->Draw("colz same"); + c2->Modified(); + c2->Update(); + + auto c23 = c.cd(4); + c23->SetGrid(1,1); + TH1F *hdiff2 = (TH1F*)(f->Get("d_sx3phi_minus_pcphi1")); + if(hdiff2) { + hdiff2->SetLineColor(kRed); + hdiff2->Draw(""); + } + TH1F *hdiff = (TH1F*)(f->Get("d_phiPC_phiQQQ_TimeCut")); + if(hdiff) { + hdiff->SetLineColor(kBlue); + hdiff->Draw("SAME"); + } + + c23->Modified(); + c23->Update(); + while(gPad->WaitPrimitive()); + + files.emplace_back(f); + std::cout <<"Test" << std::endl; + if(i==21) { + i=11; + c.Clear(); + c.Divide(2,2); + ctr=0; + } + //if(quit) break; + } + for(auto file : files) { + file->Close(); + } +} diff --git a/scratch/sx3z_vs_phiz/scan_offset_fix_witha1c1.C b/scratch/sx3z_vs_phiz/scan_offset_fix_witha1c1.C new file mode 100644 index 0000000..9916e3f --- /dev/null +++ b/scratch/sx3z_vs_phiz/scan_offset_fix_witha1c1.C @@ -0,0 +1,88 @@ +#include "testmodel.h" + +int quit=0; +void handler(int){quit=1;} + +int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5}; +void scan_offset_fix_witha1c1(){ + signal(SIGINT,handler); + TCanvas c("c1","c1",0,0,1600,800); + c.Divide(2,1); + + TF1 f1("model",model,-200,200,2); + f1.SetNpx(10000); + std::vector pars = {0.0,1.}; + f1.SetParameters(pars.data()); + f1.SetLineColor(kGreen+2); + f1.SetLineStyle(kLine); + + + + + TFile* f=NULL; + std::vector files; + int ctr=0; + for(int i=12; i<=21; i++) { + if(i==15) continue; + 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("pczfix_vs_qqqpczguess_A1C2")); + h23->SetLineColorAlpha(kOrange,0.75); + h23->Draw("box SAME"); + + TH2F *h23a = (TH2F*)(f->Get("pczfix_vs_qqqpczguess_A1C1")); + h23a->SetLineColorAlpha(kOrange,0.75); + h23a->Draw("box SAME"); + + } 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("pczfix_vs_sx3pczguess_A1C2")); + TH2F *h2a = (TH2F*)(f->Get("pczfix_vs_sx3pczguess_A1C1")); + //TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(!h2 || !h2a) continue; + h2->SetTitle(Form("case%d",i)); + //h2->Draw("colz same"); + h2->SetLineColorAlpha(colors[ctr],0.75); + h2->Draw("box same"); + h2a->SetLineColorAlpha(colors[ctr],0.75); + h2a->Draw("box same"); + //f1.Draw("same"); + } + TF1 eqline("x","x",-200,200); + eqline.Draw("SAME"); + c1->Modified(); + c1->Update(); + ctr+=1; + + + auto c2=c.cd(2); + c2->SetGrid(1,1); + + TH2F *h3 = (TH2F*)(f->Get("sx3phi_vs_pcphi1")); +// TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(!h3) continue; + h3->SetTitle(Form("case%d",i)); + h3->Draw("colz"); + eqline.Draw("SAME"); + c2->Modified(); + c2->Update(); + + while(gPad->WaitPrimitive()); + + files.emplace_back(f); + if(i==21) { + i=11; + c.Clear(); + c.Divide(2,1); + ctr=0; + } + if(quit) break; + } + for(auto file : files) { + file->Close(); + } +} diff --git a/scratch/sx3z_vs_phiz/scan_offset_witha1c1.C b/scratch/sx3z_vs_phiz/scan_offset_witha1c1.C new file mode 100644 index 0000000..636864c --- /dev/null +++ b/scratch/sx3z_vs_phiz/scan_offset_witha1c1.C @@ -0,0 +1,87 @@ +#include "testmodel.h" + +int quit=0; +void handler(int){quit=1;} + +int colors[] = {kSpring+3, kRed, kGreen+3, kBlue+3, kViolet, kOrange, kSpring-7, kAzure-5}; +void scan_offset_witha1c1(){ + signal(SIGINT,handler); + TCanvas c("c1","c1",0,0,1600,800); + c.Divide(2,1); + + TF1 f1("model",model,-200,200,2); + f1.SetNpx(10000); + std::vector pars = {0.0,1.}; + f1.SetParameters(pars.data()); + f1.SetLineColor(kGreen+2); + f1.SetLineStyle(kLine); + + + + + TFile* f=NULL; + std::vector files; + int ctr=0; + for(int i=12; i<=21; i++) { + 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")); + TH2F *h23a = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int_A1C1")); + h23->SetLineColorAlpha(kOrange,0.75); + h23->Draw("box SAME"); + h23a->SetLineColorAlpha(kOrange,0.75); + h23a->Draw("box SAME"); + + } 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")); + TH2F *h2a = (TH2F*)(f->Get("pcz_vs_sx3pczguess_A1C1")); + //TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(!h2 || !h2a) continue; + h2->SetTitle(Form("case%d",i)); + //h2->Draw("colz same"); + h2->SetLineColorAlpha(colors[ctr],0.75); + h2a->SetLineColorAlpha(colors[ctr],0.75); + h2->Draw("box same"); + h2a->Draw("box same"); + f1.Draw("same"); + } + TF1 eqline("x","x",-200,200); + eqline.Draw("SAME"); + c1->Modified(); + c1->Update(); + ctr+=1; + + + auto c2=c.cd(2); + c2->SetGrid(1,1); + + TH2F *h3 = (TH2F*)(f->Get("sx3phi_vs_pcphi1")); +// TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2")); + if(!h3) continue; + h3->SetTitle(Form("case%d",i)); + h3->Draw("colz"); + eqline.Draw("SAME"); + c2->Modified(); + c2->Update(); + + while(gPad->WaitPrimitive()); + + files.emplace_back(f); + if(i==21) { + i=11; + c.Clear(); + c.Divide(2,1); + ctr=0; + } + if(quit) break; + } + for(auto file : files) { + file->Close(); + } +} diff --git a/scratch/sx3z_vs_phiz/testmodel.h b/scratch/sx3z_vs_phiz/testmodel.h index a3deffd..1c327fe 100644 --- a/scratch/sx3z_vs_phiz/testmodel.h +++ b/scratch/sx3z_vs_phiz/testmodel.h @@ -18,7 +18,7 @@ double model_invert(double *y, double *q) { 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*3; - return result+40; + return result; } /*void testmodel() {