rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

chore: Adjust range in Analysis() function

Browse Source
This commit is contained in:
vs19g 2024-06-27 14:26:47 -04:00
parent b0093563d5
commit 862c4012eb
9 changed files with 1140 additions and 144 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Analysis.C
View File

@ -16,5 +16,5 @@ void Analysis(int start, int end) {
// Define a macro with the same name as the script // Define a macro with the same name as the script
void Analysis() { void Analysis() {
Analysis(72, 194); // Adjust the range if needed Analysis(150, 194); // Adjust the range if needed
} }

37
Analyzer.C
View File

@ -58,18 +58,17 @@ void Analyzer::Begin(TTree * /*tree*/){
hqqqVpcIndex->SetNdivisions( -612, "x"); hqqqVpcIndex->SetNdivisions( -612, "x");
hqqqVpcIndex->SetNdivisions( -12, "y"); hqqqVpcIndex->SetNdivisions( -12, "y");
hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 400, 0, 5000); hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 800, 0, 16000);
hqqqVpcE->SetNdivisions( -612, "x"); hqqqVpcE->SetNdivisions( -612, "x");
hqqqVpcE->SetNdivisions( -12, "y"); hqqqVpcE->SetNdivisions( -12, "y");
hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 400, 0, 5000); hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 800, 0, 16000);
hsx3VpcE->SetNdivisions( -612, "x"); hsx3VpcE->SetNdivisions( -612, "x");
hsx3VpcE->SetNdivisions( -12, "y"); hsx3VpcE->SetNdivisions( -12, "y");
hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600); hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000); hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000);
sx3_contr.ConstructGeo(); sx3_contr.ConstructGeo();
pw_contr.ConstructGeo(); pw_contr.ConstructGeo();
@ -283,7 +282,6 @@ Bool_t Analyzer::Process(Long64_t entry){
float aE = 0; float aE = 0;
float cE = 0; float cE = 0;
bool multi_an =false;
// if( ID[0].first < 1 ) { // if( ID[0].first < 1 ) {
// aID = pc.ch[ID[0].second]; // aID = pc.ch[ID[0].second];
// cID = pc.ch[ID[1].second]; // cID = pc.ch[ID[1].second];
@ -293,21 +291,28 @@ Bool_t Analyzer::Process(Long64_t entry){
// } // }
// printf("anode= %d, cathode = %d\n", aID, cID); // printf("anode= %d, cathode = %d\n", aID, cID);
// for( int k = 0; k < qqq.multi; k++){ for( int k = 0; k < qqq.multi; k++){
// if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){ if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){
int multi_an =0;
for(int l=0;l<E.size();l++){ for(int l=0;l<E.size();l++){
if(E[l].first<24 && E[l].first!=20 && E[l].first!=12){ if(E[l].first<24){
if(!multi_an){ multi_an++;
aE = E[l].second;
}
multi_an=true;
}
else {
cE = E[l].second + cE;
} }
} }
// }
// } if(multi_an>=1){
for(int l=0;l<E.size();l++){
if(E[l].first<24 && E[l].first!=19 && E[l].first!=12){
aE = E[l].second;
}else if(E[l].first>24){
cE = E[l].second;
}
}
}
}
}
hanVScatsum->Fill(aE,cE); hanVScatsum->Fill(aE,cE);
if( ID[0].first < 1 ) { if( ID[0].first < 1 ) {

283
Armory/#ClassPW.h# Normal file
View File

@ -0,0 +1,283 @@
#ifndef ClassPW_h
#define ClassPW_h
#include <cstdio>
#include <TMath.h>
#include <TVector3.h>
struct PWHitInfo{
std::pair<short, short> nearestWire; // anode, cathode
std::pair<double, double> nearestDist; // anode, cathode
std::pair<short, short> nextNearestWire; // anode, cathode
std::pair<double, double> nextNearestDist; // anode, cathode
void Clear(){
nearestWire.first = -1;
nearestWire.second = -1;
nearestDist.first = 999999999;
nearestDist.second = 999999999;
nextNearestWire.first = -1;
nextNearestWire.second = -1;
nextNearestDist.first = 999999999;
nextNearestDist.second = 999999999;
}
};
//!########################################################
class PW{ // proportional wire
public:
PW(){ ClearHitInfo();};
~PW(){};
PWHitInfo GetHitInfo() const {return hitInfo;}
std::pair<short, short> GetNearestID() const {return hitInfo.nearestWire;}
std::pair<double, double> GetNearestDistance() const {return hitInfo.nearestDist;}
std::pair<short, short> Get2ndNearestID() const {return hitInfo.nextNearestWire;}
std::pair<double, double> Get2ndNearestDistance() const {return hitInfo.nextNearestDist;}
TVector3 GetTrackPos() const {return trackPos;}
TVector3 GetTrackVec() const {return trackVec;}
double GetTrackTheta() const {return trackVec.Theta();}
double GetTrackPhi() const {return trackVec.Phi();}
double GetZ0();
int GetNumWire() const {return nWire;}
double GetDeltaAngle() const {return dAngle;}
double GetAnodeLength() const {return anodeLength;}
double GetCathodeLength() const {return cathodeLength;}
TVector3 GetAnodeDn(short id) const {return An[id].first;}
TVector3 GetAnodeUp(short id) const {return An[id].second;}
TVector3 GetCathodeDn(short id) const {return Ca[id].first;}
TVector3 GetCathodeUp(short id) const {return Ca[id].second;}
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();}
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 FindWireID(TVector3 pos, TVector3 direction, bool verbose = false);
void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false);
void CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false);
void Print(){
printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first,
hitInfo.nearestDist.first,
hitInfo.nearestWire.second,
hitInfo.nearestDist.second);
printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire.first,
hitInfo.nextNearestDist.first,
hitInfo.nextNearestWire.second,
hitInfo.nextNearestDist.second);
}
private:
PWHitInfo hitInfo;
TVector3 trackPos;
TVector3 trackVec;
const int nWire = 24;
const int wireShift = 3;
const float zLen = 380; //mm
const float radiusA = 37;
const float radiusC = 43;
double dAngle;
double anodeLength;
double cathodeLength;
std::vector<std::pair<TVector3,TVector3>> An; // the anode wire position vector in space
std::vector<std::pair<TVector3,TVector3>> Ca; // the cathode wire position vector in space
double Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2){
TVector3 na = a1 - a2;
TVector3 nb = b1 - b2;
TVector3 nd = (na.Cross(nb)).Unit();
return TMath::Abs(nd.Dot(a1-b2));
}
};
inline void PW::ClearHitInfo(){
hitInfo.Clear();
}
inline void PW::ConstructGeo(){
An.clear();
Ca.clear();
std::pair<TVector3, TVector3> p1; // anode
std::pair<TVector3, TVector3> q1; // cathode
//anode and cathode start at pos-Y axis and count in right-Hand
//anode wire shift is right-hand.
//cathode wire shift is left-hand.
for(int i = 0; i < nWire; i++ ){
// Anode rotate right-hand
p1.first.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
radiusA * TMath::Sin( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
zLen/2);
p1.second.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
radiusA * TMath::Sin( TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
-zLen/2);
An.push_back(p1);
// Cathod rotate left-hand
q1.first.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
radiusC * TMath::Sin( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
zLen/2);
q1.second.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()),
radiusC * TMath::Sin( TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()),
-zLen/2);
Ca.push_back(q1);
}
dAngle = wireShift * TMath::TwoPi() / nWire;
anodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusA * TMath::Sin(dAngle/2),2) );
cathodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusC * TMath::Sin(dAngle/2),2) );
}
inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose ){
hitInfo.Clear();
double phi = direction.Phi();
for( int i = 0; i < nWire; i++){
double disA = 99999999;
double phiS = An[i].first.Phi() - TMath::PiOver4();
double phiL = An[i].second.Phi() + TMath::PiOver4();
// printf("A%2d: %f %f | %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg(), phi * TMath::RadToDeg());
if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi();
if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi();
if( phiS < phi && phi < phiL) {
disA = Distance( pos, pos + direction, An[i].first, An[i].second);
if( disA < hitInfo.nearestDist.first ){
hitInfo.nearestDist.first = disA;
hitInfo.nearestWire.first = i;
}
}
double disC = 99999999;
phiS = Ca[i].second.Phi()- TMath::PiOver4();
phiL = Ca[i].first.Phi() + TMath::PiOver4();
// printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi();
if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi();
if(phiS < phi && phi < phiL) {
disC = Distance( pos, pos + direction, Ca[i].first, Ca[i].second);
if( disC < hitInfo.nearestDist.second ){
hitInfo.nearestDist.second = disC;
hitInfo.nearestWire.second = i;
}
}
if(verbose) printf(" %2d | %8.2f, %8.2f\n", i, disA, disC);
}
//==== find the 2nd nearest wire
short anode1 = hitInfo.nearestWire.first;
short aaa1 = anode1 - 1; if( aaa1 < 0 ) aaa1 += nWire;
short aaa2 = (anode1 + 1) % nWire;
double haha1 = Distance( pos, pos + direction, An[aaa1].first, An[aaa1].second);
double haha2 = Distance( pos, pos + direction, An[aaa2].first, An[aaa2].second);
if( haha1 < haha2){
hitInfo.nextNearestWire.first = aaa1;
hitInfo.nextNearestDist.first = haha1;
}else{
hitInfo.nextNearestWire.first = aaa2;
hitInfo.nextNearestDist.first = haha2;
}
short cathode1 = hitInfo.nearestWire.second;
short ccc1 = cathode1 - 1; if( ccc1 < 0 ) ccc1 += nWire;
short ccc2 = (cathode1 + 1) % nWire;
haha1 = Distance( pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second);
haha2 = Distance( pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second);
if( haha1 < haha2){
hitInfo.nextNearestWire.second = ccc1;
hitInfo.nextNearestDist.second = haha1;
}else{
hitInfo.nextNearestWire.second = ccc2;
hitInfo.nextNearestDist.second = haha2;
}
if( verbose ) Print();
}
inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose){
trackPos = sx3Pos;
TVector3 n1 = (An[anodeID].first - An[anodeID].second).Cross((sx3Pos - An[anodeID].second)).Unit();
TVector3 n2 = (Ca[cathodeID].first - Ca[cathodeID].second).Cross((sx3Pos - Ca[cathodeID].second)).Unit();
// if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit();
if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg());
}
inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose){
trackPos = sx3Pos;
double p1 = TMath::Abs(hitInfo.nearestDist.first + gRandom->Gaus(0, sigmaA));
double p2 = TMath::Abs(hitInfo.nextNearestDist.first + gRandom->Gaus(0, sigmaA));
double fracA = p1 / (p1 + p2);
short anodeID1 = hitInfo.nearestWire.first;
short anodeID2 = hitInfo.nextNearestWire.first;
TVector3 shiftA1 = (An[anodeID2].first - An[anodeID1].first) * fracA;
TVector3 shiftA2 = (An[anodeID2].second - An[anodeID1].second) * fracA;
double q1 = TMath::Abs(hitInfo.nearestDist.second + gRandom->Gaus(0, sigmaC));
double q2 = TMath::Abs(hitInfo.nextNearestDist.second + gRandom->Gaus(0, sigmaC));
double fracC = q1 / (q1 + q2);
short cathodeID1 = hitInfo.nearestWire.second;
short cathodeID2 = hitInfo.nextNearestWire.second;
TVector3 shiftC1 = (Ca[cathodeID2].first - Ca[cathodeID1].first) * fracC;
TVector3 shiftC2 = (Ca[cathodeID2].second - Ca[cathodeID1].second) * fracC;
TVector3 a1 = An[anodeID1].first + shiftA1;
TVector3 a2 = An[anodeID1].second + shiftA2;
TVector3 c1 = Ca[cathodeID1].first + shiftC1;
TVector3 c2 = Ca[cathodeID1].second + shiftC2;
TVector3 n1 = (a1 - a2).Cross((sx3Pos - a2)).Unit();
TVector3 n2 = (c1 - c2).Cross((sx3Pos - c2)).Unit();
// if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit();
if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg());
}
inline double PW::GetZ0(){
double x = trackPos.X();
double y = trackPos.Y();
double rho = TMath::Sqrt(x*x + y*y);
double theta = trackVec.Theta();
return trackPos.Z() - rho / TMath::Tan(theta);
}
#endif

20
Armory/ClassPW.h
View File

@ -4,6 +4,7 @@
#include <cstdio> #include <cstdio>
#include <TMath.h> #include <TMath.h>
#include <TVector3.h> #include <TVector3.h>
#include <TRandom.h>
struct PWHitInfo{ struct PWHitInfo{
std::pair<short, short> nearestWire; // anode, cathode std::pair<short, short> nearestWire; // anode, cathode
@ -65,8 +66,6 @@ public:
void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false); void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false);
void CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false); void CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false);
double CircularMean(std::vector<std::pair<int, double>> wireList);
void Print(){ void Print(){
printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first, printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first,
hitInfo.nearestDist.first, hitInfo.nearestDist.first,
@ -282,21 +281,4 @@ inline double PW::GetZ0(){
} }
inline double PW::CircularMean(std::vector<std::pair<int, double>> wireList){
//use unit vector, wireID start from Zero
double xCom = 0, yCom = 0;
for( size_t i = 0; i < wireList.size() ; i++){
xCom += TMath::Cos(TMath::TwoPi() * wireList[i].first / nWire) * wireList[i].second;
yCom += TMath::Sin(TMath::TwoPi() * wireList[i].first / nWire) * wireList[i].second;
}
//calculate the angle of the summed unit vectors
double angle = TMath::ATan2(yCom, xCom);
if( angle < 0 ) angle += TMath::TwoPi(); // convert the angle from 0 to 2 pi
return angle/ TMath::TwoPi() * nWire;
}
#endif #endif

158
TrackRecon.C
View File

@ -1,6 +1,6 @@
#define Analyzer_cxx #define TrackRecon_cxx
#include "Analyzer.h" #include "TrackRecon.h"
#include <TH2.h> #include <TH2.h>
#include <TStyle.h> #include <TStyle.h>
#include <TCanvas.h> #include <TCanvas.h>
@ -14,20 +14,6 @@
#include "TVector3.h" #include "TVector3.h"
TH2F * hsx3IndexVE;
TH2F * hqqqIndexVE;
TH2F * hpcIndexVE;
TH2F * hsx3Coin;
TH2F * hqqqCoin;
TH2F * hpcCoin;
TH2F * hqqqPolar;
TH2F * hsx3VpcIndex;
TH2F * hqqqVpcIndex;
TH2F * hqqqVpcE;
TH2F * hsx3VpcE;
TH2F * hanVScatsum;
int padID = 0; int padID = 0;
SX3 sx3_contr; SX3 sx3_contr;
@ -37,7 +23,7 @@ bool HitNonZero;
TH1F * hZProj; TH1F * hZProj;
void Analyzer::Begin(TTree * /*tree*/){ void TrackRecon::Begin(TTree * /*tree*/){
TString option = GetOption(); TString option = GetOption();
hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600); hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
@ -47,17 +33,14 @@ void Analyzer::Begin(TTree * /*tree*/){
} }
Bool_t TrackRecon::Process(Long64_t entry){
Bool_t Analyzer::Process(Long64_t entry){
// if ( entry > 100 ) return kTRUE; // if ( entry > 100 ) return kTRUE;
hitPos.Clear(); hitPos.Clear();
HitNonZero = false; HitNonZero = false;
// if( entry > 1) return kTRUE; if( entry > 1) return kTRUE;
// printf("################### ev : %llu \n", entry); // printf("################### ev : %llu \n", entry);
b_sx3Multi->GetEntry(entry); b_sx3Multi->GetEntry(entry);
@ -133,13 +116,6 @@ Bool_t Analyzer::Process(Long64_t entry){
}else{ }else{
sx3ChBk = sx3.ch[index]; sx3ChBk = sx3.ch[index];
} }
for( int j = 0; j < pc.multi; j++){
// hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
if( sx3.ch[index] > 8 ){
hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
// hpcIndexVE->Fill( pc.index[i], pc.e[i] );
}
}
} }
sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn); sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
@ -170,7 +146,7 @@ Bool_t Analyzer::Process(Long64_t entry){
double theta = -TMath::Pi()/2 + 2*TMath::Pi()/16/4.*(qqq.id[i]*16 + chWedge +0.5); double theta = -TMath::Pi()/2 + 2*TMath::Pi()/16/4.*(qqq.id[i]*16 + chWedge +0.5);
double rho = 10.+40./16.*(chRing+0.5); double rho = 10.+40./16.*(chRing+0.5);
// if(qqq.e[i]>50){ // if(qqq.e[i]>50){
hqqqPolar->Fill( theta, rho); // hqqqPolar->Fill( theta, rho);
// } // }
// qqq.used[i] = true; // qqq.used[i] = true;
// qqq.used[j] = true; // qqq.used[j] = true;
@ -187,7 +163,6 @@ Bool_t Analyzer::Process(Long64_t entry){
} }
// //======================= PC // //======================= PC
PCHit_1An hitInfo;
ID.clear(); ID.clear();
int counter=0; int counter=0;
@ -197,33 +172,51 @@ Bool_t Analyzer::Process(Long64_t entry){
if( E.size()==3 ){ if( E.size()==3 ){
float aE = 0; float aE = 0;
float cE = 0; float cE = 0;
bool multi_an =false; int multi_an =0;
for(int l=0;l<E.size();l++){ for(int l=0;l<E.size();l++){
if(E[l].first<24 && E[l].first!=20 && E[l].first!=12){ if(E[l].first<24 && E[l].first!=20 && E[l].first!=12){
if(!multi_an){ multi_an++;
aE = E[l].second;
}
multi_an=true;
}
else {
cE = E[l].second + cE;
}
} }
// printf("anode= %d, cathode = %d\n", aID, cID);
// }
if( ID[0].first < 1 ) {
aID = pc.ch[ID[0].second];
cID = pc.ch[ID[1].second];
}else{
cID = pc.ch[ID[0].second];
aID = pc.ch[ID[1].second];
} }
hanVScatsum->Fill(aE,cE); if(multi_an==1){
for(int l=0;l<E.size();l++){
if(E[l].first<24 && E[l].first!=20 && E[l].first!=12){
aE = E[l].second;
}else if(E[l].first>24){
cE = E[l].second;
}
}
}
//using CalTrack3 to get the track position and direction
// hanVScatsum->Fill(aE,cE);
if( HitNonZero){ if( HitNonZero){
pw_contr.CalTrack3( hitPos, hitinfo, cID); if (ID.size() == 3) {
hZProj->Fill(pw_contr.GetZ0()); int aID = -1;
int cID1 = -1;
int cID2 = -1;
for (int i = 0; i < ID.size(); i++) {
if (pc.ch[ID[i].second] < 24 && pc.ch[ID[i].second] != 20 && pc.ch[ID[i].second] != 12) {
aID = pc.ch[ID[i].second];
} else if (pc.ch[ID[i].second] > 24) {
if (cID1 == -1) {
cID1 = pc.ch[ID[i].second];
} else {
cID2 = pc.ch[ID[i].second];
}
}
}
if (aID != -1 && cID1 != -1 && cID2 != -1) {
pw_contr.CalTrack3(hitPos, aID, cID1, cID2);
pw_contr.Print();
printf("###################\n");
hZProj->Fill(pw_contr.GetZ0());
}
}
} }
// } // }
@ -240,64 +233,13 @@ Bool_t Analyzer::Process(Long64_t entry){
return kTRUE; return kTRUE;
} }
void Analyzer::Terminate(){ void TrackRecon::Terminate(){
gStyle->SetOptStat("neiou"); gStyle->SetOptStat("neiou");
TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000); TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 200, 200);
canvas->Divide(3,3); padID=1;
//hsx3VpcIndex->Draw("colz");
//=============================================== pad-1
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hsx3IndexVE->Draw("colz");
//=============================================== pad-2
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hqqqIndexVE->Draw("colz");
//=============================================== pad-3
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hpcIndexVE->Draw("colz");
//=============================================== pad-4
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hsx3Coin->Draw("colz");
//=============================================== pad-5
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hqqqCoin->Draw("colz");
//=============================================== pad-6
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hpcCoin->Draw("colz");
//=============================================== pad-7
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hsx3VpcIndex ->Draw("colz");
// hsx3VpcE->Draw("colz") ;
//=============================================== pad-8
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hqqqVpcIndex ->Draw("colz");
// hqqqVpcE ->Draw("colz");
//=============================================== pad-9
padID ++;
// canvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
// hqqqPolar->Draw("same colz pol");
canvas->cd(padID); canvas->cd(padID)->SetGrid(1); canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hZProj->Draw(); hZProj->Draw();
// hanVScatsum->Draw("colz");
} }

114
TrackRecon.h Normal file
View File

@ -0,0 +1,114 @@
#ifndef TrackRecon_h
#define TrackRecon_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include "Armory/ClassDet.h"
class TrackRecon : public TSelector {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
// Fixed size dimensions of array or collections stored in the TTree if any.
// Declaration of leaf types
Det sx3;
Det qqq;
Det pc ;
ULong64_t evID;
UInt_t run;
// List of branches
TBranch *b_eventID; //!
TBranch *b_run; //!
TBranch *b_sx3Multi; //!
TBranch *b_sx3ID; //!
TBranch *b_sx3Ch; //!
TBranch *b_sx3E; //!
TBranch *b_sx3T; //!
TBranch *b_qqqMulti; //!
TBranch *b_qqqID; //!
TBranch *b_qqqCh; //!
TBranch *b_qqqE; //!
TBranch *b_qqqT; //!
TBranch *b_pcMulti; //!
TBranch *b_pcID; //!
TBranch *b_pcCh; //!
TBranch *b_pcE; //!
TBranch *b_pcT; //!
TrackRecon(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~TrackRecon() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree);
virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
virtual void SetOption(const char *option) { fOption = option; }
virtual void SetObject(TObject *obj) { fObject = obj; }
virtual void SetInputList(TList *input) { fInput = input; }
virtual TList *GetOutputList() const { return fOutput; }
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(TrackRecon,0);
};
#endif
#ifdef TrackRecon_cxx
void TrackRecon::Init(TTree *tree){
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_eventID);
fChain->SetBranchAddress("run", &run, &b_run);
sx3.SetDetDimension(24,12);
qqq.SetDetDimension(4,32);
pc.SetDetDimension(2,24);
fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID);
fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch);
fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E);
fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T);
fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti);
fChain->SetBranchAddress("qqqID", &qqq.id, &b_qqqID);
fChain->SetBranchAddress("qqqCh", &qqq.ch, &b_qqqCh);
fChain->SetBranchAddress("qqqE", &qqq.e, &b_qqqE);
fChain->SetBranchAddress("qqqT", &qqq.t, &b_qqqT);
fChain->SetBranchAddress("pcMulti", &pc.multi, &b_pcMulti);
fChain->SetBranchAddress("pcID", &pc.id, &b_pcID);
fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh);
fChain->SetBranchAddress("pcE", &pc.e, &b_pcE);
fChain->SetBranchAddress("pcT", &pc.t, &b_pcT);
}
Bool_t TrackRecon::Notify(){
return kTRUE;
}
void TrackRecon::SlaveBegin(TTree * /*tree*/){
TString option = GetOption();
}
void TrackRecon::SlaveTerminate(){
}
#endif // #ifdef TrackRecon_cxx

15
archivist Normal file
View File

@ -0,0 +1,15 @@
#!/bin/bash
RUNNO=$1
BINARYDIR=/media/nvmeData/ANASEN_ppprimeMar20/
ARCHIVE=/media/nvmeData/ANASEN_ppprimeMar20/Archive/run_$RUNNO.tar.gz
echo "Running archivist for binary data in $BINARYDIR to archive $ARCHIVE..."
cd $BINARYDIR
tar -cvzf $ARCHIVE ./*Run_$RUNNO*.fsu
cd -
echo "Complete."

541
gainmatch.C Normal file
View File

@ -0,0 +1,541 @@
#define gainmatch_cxx
#include "gainmatch.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TCutG.h>
#include <utility>
#include <algorithm>
#include "Armory/ClassSX3.h"
#include "Armory/ClassPW.h"
#include "TVector3.h"
TH2F * hsx3IndexVE;
TH2F * hqqqIndexVE;
TH2F * hqqqIndexVE_cut;
TH2F * hpcIndexVE;
TH2F * hsx3Coin;
TH2F * hqqqCoin;
TH2F * hpcCoin;
TH2F * hpcCoin_cut;
TH2F * hGoodQQQ;
TH2F * hGoodQQQRingVWedge;
TH2F * hqqqPolar;
TH2F * hsx3VpcIndex;
TH2F * hqqqVpcIndex;
TH2F * hqqqVpcIndex_cut;
TH2F * hqqqVpcE;
TH2F * hqqqVpcE_cut;
TH2F * hqqqVpcE_cut1;
TH2F * hqqqVpcE_cut2;
TH2F * hqqqVpcE_cutCoinc;
TH2F * hsx3VpcE;
TH2F * hanVScatsum;
TH2F * hanVScatsum_cut;
TH2F * hanVScatsum_cut1;
TH2F * hanVScatsum_cut2;
TH2F * hsx3Vsx3;
TH2F * hsx3uVsx3d_01;
TH2F * hsx3uVsx3d_23;
TH2F * hsx3uVsx3d_45;
TH2F * hsx3uVsx3d_67;
TH2F * hVCID;
TH1F *hsx3bk_9_shifted ;
TH1F *hsx3bk_10_shifted ;
TH1F *hsx3bk_11_shifted ;
int padID = 0;
TCutG *Coinc_cut_set1;
//TCutG *crap_cut;
TCutG *AnCathCoinc_cut;
TCutG *AnCathCoinc_cut1;
TCutG *AnCathCoinc_cut2;
SX3 sx3_contr;
PW pw_contr;
TVector3 hitPos;
bool HitNonZero;
bool inCut;
bool inCut1;
bool inCut2;
bool inCutCoinc;
TH1F *hZd_01_1;
TH1F *hZd_01_2;
TH1F *hZd_01_3;
TH1F *hZd_01_4;
TH1F * hZProj;
TH1F * hsx3bk_11;
TH1F * hsx3bk_10;
TH1F * hsx3bk_9;
TH1F * hsx3bk_8;
void gainmatch::Begin(TTree * /*tree*/){
TString option = GetOption();
hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24*12, 0, 24*12, 400, 0, 5000); hsx3IndexVE->SetNdivisions( -612, "x");
hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
hqqqIndexVE_cut = new TH2F("hqqqIndexVE_cut", "QQQ index vs Energy gated; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2*24, 0, 2*24, 400, 0, 4000); hpcIndexVE->SetNdivisions( -1204, "x");
hGoodQQQ = new TH2F("hGoodQQQ", "number of good QQQ vs QQQ id", 10, 0, 10, 4, 0, 4);
hGoodQQQRingVWedge = new TH2F("hGoodQQQRingVWedge", "Ring index, Wedge index", 16*4, 0, 16*4, 16*4, 0, 16*4);
hZd_01_1 =new TH1F("hZd_01_1", "Z position", 100, -1, 1);
hZd_01_2 =new TH1F("hZd_01_2", "Z position", 100, -1, 1);
hZd_01_3 =new TH1F("hZd_01_3", "Z position", 100, -1, 1);
hZd_01_4 =new TH1F("hZd_01_4", "Z position", 100, -1, 1);
hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24*12, 0, 24*12, 24*12, 0, 24*12);
hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4*2*16, 0, 4*2*16, 4*2*16, 0, 4*2*16);
hpcCoin = new TH2F("hpcCoin", "PC Coincident", 2*24, 0, 2*24, 2*24, 0, 2*24);
hpcCoin_cut = new TH2F("hpcCoin_cut", "PC Coincident gated", 2*24, 0, 2*24, 2*24, 0, 2*24);
hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16*4, -TMath::Pi(), TMath::Pi(),16, 10, 50);
hsx3VpcIndex = new TH2F("hsx3Vpcindex", "sx3 vs pc; sx3 index; pc index", 24*12, 0, 24*12, 48, 0, 48);
hsx3Vsx3 = new TH2F("hsx3Vsx3", "sx3 vs sx3; sx3 E; sx3 E", 8000, 0, 16000, 8000, 0, 16000);
hsx3uVsx3d_01 = new TH2F("hsx3uVsx3d_01", "sx3u vs sx3d; sx3u E; sx3d E", 100, 0, 1, 100, 0, 1);
hsx3uVsx3d_23 = new TH2F("hsx3uVsx3d_23", "sx3u vs sx3d; sx3u E; sx3d E", 100, 0, 1, 100, 0, 1);
hsx3uVsx3d_45 = new TH2F("hsx3uVsx3d_45", "sx3u vs sx3d; sx3u E; sx3d E", 1000, 0, 5000, 1000, 0, 5000);
hsx3uVsx3d_67 = new TH2F("hsx3uVsx3d_67", "sx3u vs sx3d; sx3u E; sx3d E", 1000, 0, 5000, 1000, 0, 5000);
hsx3VpcIndex->SetNdivisions( -612, "x");
hsx3VpcIndex->SetNdivisions( -12, "y");
hqqqVpcIndex = new TH2F("hqqqVpcindex", "qqq vs pc; qqq index; pc index", 4*2*16, 0, 4*2*16, 48, 0, 48);
hqqqVpcIndex->SetNdivisions( -612, "x");
hqqqVpcIndex->SetNdivisions( -12, "y");
hqqqVpcIndex_cut = new TH2F("hqqqVpcindex_cut", "qqq vs pc gated; qqq index; pc index", 4*2*16, 0, 4*2*16, 48, 0, 48);
hqqqVpcIndex_cut->SetNdivisions( -612, "x");
hqqqVpcIndex_cut->SetNdivisions( -12, "y");
hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
hqqqVpcE->SetNdivisions( -612, "x");
hqqqVpcE->SetNdivisions( -12, "y");
hqqqVpcE_cut = new TH2F("hqqqVpcEnergy_cut", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
hqqqVpcE_cut->SetNdivisions( -612, "x");
hqqqVpcE_cut->SetNdivisions( -12, "y");
hqqqVpcE_cut1 = new TH2F("hqqqVpcEnergy_cut1", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
hqqqVpcE_cut1->SetNdivisions( -612, "x");
hqqqVpcE_cut1->SetNdivisions( -12, "y");
hqqqVpcE_cut2 = new TH2F("hqqqVpcEnergy_cut2", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
hqqqVpcE_cut2->SetNdivisions( -612, "x");
hqqqVpcE_cut2->SetNdivisions( -12, "y");
hqqqVpcE_cutCoinc = new TH2F("hqqqVpcEnergy_cutCoinc", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
hqqqVpcE_cutCoinc->SetNdivisions( -612, "x");
hqqqVpcE_cutCoinc->SetNdivisions( -12, "y");
hsx3bk_8=new TH1F("hsx3bk_8", "hsx3bk_8",1000, 0,5000);
hsx3bk_9=new TH1F("hsx3bk_9", "hsx3bk_9",1000, 0,5000);
hsx3bk_10=new TH1F("hsx3bk_10", "hsx3bk_10",1000, 0,5000);
hsx3bk_11=new TH1F("hsx3bk_11", "hsx3bk_11",1000, 0,5000);
hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 400, 0, 5000);
hsx3VpcE->SetNdivisions( -612, "x");
hsx3VpcE->SetNdivisions( -12, "y");
hsx3bk_9_shifted = new TH1F("hsx3bk_9_shifted", "hsx3bk_9",1000, 0,5000);
hsx3bk_10_shifted = new TH1F("hsx3bk_10_shifted", "hsx3bk_9",1000, 0,5000);
hsx3bk_11_shifted = new TH1F("hsx3bk_11_shifted", "hsx3bk_9",1000, 0,5000);
hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 8000,0 , 16000, 8000, 0 , 16000);
hanVScatsum_cut = new TH2F("hanVScatsum_cut", "Anode vs Cathode Sum gated; Anode E; Cathode E", 1600,0 , 16000, 1600, 0 , 16000);
hanVScatsum_cut1 = new TH2F("hanVScatsum_cut1", "Anode vs Cathode Sum gated; Anode E; Cathode E", 1600,0 , 16000, 1600, 0 , 16000);
hanVScatsum_cut2 = new TH2F("hanVScatsum_cut2", "Anode vs Cathode Sum gated; Anode E; Cathode E", 1600,0 , 16000, 1600, 0 , 16000);
hVCID = new TH2F("hVCID", "Virtual Cathod ID vs total Cath. Energy", 200, 0, 24, 200, 0, 10000);
sx3_contr.ConstructGeo();
pw_contr.ConstructGeo();
TFile *f3 = new TFile("Coinc_cut_set1.root");
//TFile *f4 = new TFile("crap_cut.root");
TFile *f = new TFile("AnCathCoinc_cut.root");
TFile *f1 = new TFile("AnCathCoinc_cut1.root");
TFile *f2 = new TFile("AnCathCoinc_cut2.root");
Coinc_cut_set1 = (TCutG*)f3->Get("Coinc_cut_set1");
//crap_cut = (TCutG*)f4->Get("crap_cut");
AnCathCoinc_cut = (TCutG*)f->Get("AnCathCoinc_cut");
AnCathCoinc_cut1 = (TCutG*)f1->Get("AnCathCoinc_cut1");
AnCathCoinc_cut2 = (TCutG*)f2->Get("AnCathCoinc_cut2");
}
Bool_t gainmatch::Process(Long64_t entry){
// if ( entry > 100 ) return kTRUE;
hitPos.Clear();
HitNonZero = false;
inCut = false;
// if( entry > 1) return kTRUE;
// printf("################### ev : %llu \n", entry);
b_sx3Multi->GetEntry(entry);
b_sx3ID->GetEntry(entry);
b_sx3Ch->GetEntry(entry);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
b_qqqMulti->GetEntry(entry);
b_qqqID->GetEntry(entry);
b_qqqCh->GetEntry(entry);
b_qqqE->GetEntry(entry);
b_qqqT->GetEntry(entry);
b_pcMulti->GetEntry(entry);
b_pcID->GetEntry(entry);
b_pcCh->GetEntry(entry);
b_pcE->GetEntry(entry);
b_pcT->GetEntry(entry);
sx3.CalIndex();
qqq.CalIndex();
pc.CalIndex();
// sx3.Print();
//########################################################### Raw data
// //======================= SX3
std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
for (int i = 0; i < sx3.multi; i++) {
ID.push_back(std::pair<int, int>(sx3.id[i], i));
hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
for (int j = i + 1; j < sx3.multi; j++) {
hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
}
for (int j = 0; j < pc.multi; j++) {
hsx3VpcIndex->Fill(sx3.index[i], pc.index[j]);
}
}
if (ID.size() > 0) {
std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> &a, const std::pair<int, int> &b) {
return a.first < b.first;
});
std::vector<std::pair<int, int>> sx3ID;
sx3ID.push_back(ID[0]);
bool found = false;
for (size_t i = 1; i < ID.size(); i++) {
if (ID[i].first == sx3ID.back().first) {
sx3ID.push_back(ID[i]);
if (sx3ID.size() >= 3) {
found = true;
}
} else {
if (!found) {
sx3ID.clear();
sx3ID.push_back(ID[i]);
}
}
}
if (found) {
int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1;
float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0;
for (size_t i = 0; i < sx3ID.size(); i++) {
int index = sx3ID[i].second;
if (sx3.ch[index] < 8) {
if (sx3.ch[index] % 2 == 0) {
sx3ChDn = sx3.ch[index];
sx3EDn = sx3.e[index];
} else {
sx3ChUp = sx3.ch[index];
sx3EUp = sx3.e[index];
}
} else {
sx3ChBk = sx3.ch[index];
sx3EBk = sx3.e[index];
}
int ch = sx3.ch[index];
float energy = sx3.e[index];
if (sx3ID[0].first == 9) {
float peak8 = 0.0;
float peak9 = 0.0;
int peak10 = 0.0;
float peak11 = 0.0;
float shift9 =0.0;
float shift10 =0.0;
float shift11 =0.0;
int minBin_8 = hsx3bk_8->FindBin(1);
int maxBin_8 = hsx3bk_8->FindBin(5000);
int maxRangeBinContent_8 = -1;
double maxBinCenter_8 = 0.0;
int minBin_9 = hsx3bk_9->FindBin(1);
int maxBin_9 = hsx3bk_9->FindBin(5000);
int maxRangeBinContent_9 = -1;
double maxBinCenter_9 = 0.0;
int minBin_10 = hsx3bk_10->FindBin(1);
int maxBin_10 = hsx3bk_10->FindBin(5000);
int maxRangeBinContent_10 = -1;
double maxBinCenter_10 = 0.0;
int minBin_11 = hsx3bk_11->FindBin(1);
int maxBin_11 = hsx3bk_11->FindBin(5000);
int maxRangeBinContent_11 = -1;
double maxBinCenter_11 = 0.0;
if (sx3ChBk == 8) {
hsx3bk_8->Fill(sx3EBk);
for (int bin = minBin_8; bin <= maxBin_8; ++bin) {
if (hsx3bk_8->GetBinContent(bin) > maxRangeBinContent_8) {
maxRangeBinContent_8 = hsx3bk_8->GetBinContent(bin);
maxBinCenter_8 = hsx3bk_8->GetBinCenter(bin);
}
}
//peak8 = hsx3bk_8->GetMaximumBin();
//peak8 = hsx3bk_8->GetMaximumBin();
//printf("peak8: %f\n", maxBinCenter_8);
}
//printf("peak8_mm: %f\n", maxBinCenter);
else if (sx3ChBk == 9) {
hsx3bk_9->Fill(sx3EBk);
for (int bin = minBin_9; bin <= maxBin_9; ++bin) {
if (hsx3bk_9->GetBinContent(bin) > maxRangeBinContent_9) {
maxRangeBinContent_9 = hsx3bk_9->GetBinContent(bin);
maxBinCenter_9 = hsx3bk_9->GetBinCenter(bin);
}
}
//peak8 = hsx3bk_8->GetMaximumBin();
//peak8 = hsx3bk_8->GetMaximumBin();
//printf("peak9: %f\n", maxBinCenter_9);
//hsx3bk_9_shifted->Fill(sx3EBk*0.76);
peak9 = 2097.5/maxBinCenter_9;
//printf("peak9_shift: %f\n", peak9);
hsx3bk_9_shifted->Fill(sx3EBk*(2097.5/maxBinCenter_9));
//printf("peak9 %d\n", peak9);
}
else if(sx3ChBk == 10) {
hsx3bk_10->Fill(sx3EBk);
for (int bin = minBin_10; bin <= maxBin_10; ++bin) {
if (hsx3bk_10->GetBinContent(bin) > maxRangeBinContent_10) {
maxRangeBinContent_10 = hsx3bk_10->GetBinContent(bin);
maxBinCenter_10 = hsx3bk_10->GetBinCenter(bin);
}
}
//peak8 = hsx3bk_8->GetMaximumBin();
//peak8 = hsx3bk_8->GetMaximumBin();
//printf("peak10: %f\n", maxBinCenter_10);
//hsx3bk_9_shifted->Fill(sx3EBk*0.76);
peak10= 2097.5/maxBinCenter_10;
//printf("peak10_shift: %f\n", 1787.5/maxBinCenter_10);
hsx3bk_10_shifted->Fill(sx3EBk*(2097.5/maxBinCenter_10));
//printf("peak9 %d\n", peak9);
}
//peak10 = hsx3bk_10->GetMaximumBin();
// printf("peak10 %d\n" ,peak10);
else if(sx3ChBk == 11) {
hsx3bk_11->Fill(sx3EBk);
for (int bin = minBin_11; bin <= maxBin_11; ++bin) {
if (hsx3bk_11->GetBinContent(bin) > maxRangeBinContent_11) {
maxRangeBinContent_11 = hsx3bk_11->GetBinContent(bin);
maxBinCenter_11 = hsx3bk_11->GetBinCenter(bin);
}
}
//peak8 = hsx3bk_8->GetMaximumBin();
//peak8 = hsx3bk_8->GetMaximumBin();
//printf("peak9: %f\n", maxBinCenter_11);
//hsx3bk_9_shifted->Fill(sx3EBk*0.76);
peak11 = 2097.5/maxBinCenter_11;
//printf("peak11_shift: %f\n", peak11);
hsx3bk_11_shifted->Fill(sx3EBk*(2097.5/maxBinCenter_11));
//printf("peak9 %d\n", peak9);
}
float sx3EBk_shifted = 0.0;
float sx3E_u_matched_01 = 0.0;
float sx3E_d_matched_01 = 0.0;
float sx3E_fb_matched_01 = 0.0;
float sx3E_fbu_matched_01 = 0.0;
float sx3E_fbd_matched_01 = 0.0;
float diff =0.0;
float ratio = 0.0;
float coeff = 0.0;
if (sx3ChBk == 9) {
sx3EBk_shifted = (sx3EBk *(2097.5/maxBinCenter_9));
} else if (sx3ChBk == 10) {
sx3EBk_shifted = (sx3EBk * (2097.5/maxBinCenter_10));
} else if (sx3ChBk == 11) {
sx3EBk_shifted = (sx3EBk * (2097.5/maxBinCenter_11)) ;
} else {
sx3EBk_shifted = sx3EBk; // Use unshifted value for sx3ChBk == 8
}
if ((sx3ChUp == 1 && sx3ChDn == 0)) {
sx3E_u_matched_01= (sx3EUp-0.898729)/0.836243;
//sx3E_u_matched_01= (0.836243*sx3EDn)+0.898729;
sx3E_d_matched_01= (sx3EDn-0.898729)/0.836243;
sx3E_fb_matched_01=(sx3EBk_shifted+9.2423)/0.924773 ;
sx3E_fbu_matched_01=(sx3E_u_matched_01+9.2423)/0.924773 ;
sx3E_fbd_matched_01=(sx3E_d_matched_01+9.2423)/0.924773 ;
diff = sx3E_fb_matched_01 - (sx3EUp+sx3E_fbd_matched_01);
ratio = sx3EUp/sx3E_fbd_matched_01;
coeff = ((sx3EUp+diff) - (sx3E_fbd_matched_01*ratio))/(diff*(1+ratio));
}
//TH2F *hsx3uVsx3d_01 = nullptr;
if (sx3ChBk >=8) {
//if (sx3ChBk == 9) {
if ((sx3ChUp == 1 && sx3ChDn == 0)) {
if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
if (sx3EBk_shifted > 50 && sx3EUp > 50 && sx3EDn > 50) {
printf("sx3EUp: %f, sx3EDn: %f, sx3E_u_matched_01: %f,sx3E_d_matched_01: %f\n", sx3EUp, sx3EDn, sx3E_u_matched_01,sx3E_d_matched_01);
//printf("Filling hsx3uVsx3d_01_shifted: %f\n", sx3EBk_ud_matched_01 / sx3EBk_shifted);
// hsx3uVsx3d_01->Fill(sx3E_u_matched_01 / sx3EBk_shifted, sx3E_d_matched_01 / sx3EBk_shifted);
hsx3uVsx3d_01->Fill(sx3EUp / sx3EBk_shifted, sx3E_d_matched_01 / sx3EBk_shifted);
hsx3uVsx3d_23->Fill(sx3EUp / sx3EBk_shifted, sx3EDn/ sx3EBk_shifted);
}
}
//}
}
else if ((sx3ChUp == 3 && sx3ChDn == 2)) {
if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
if (sx3EBk_shifted != 0 && sx3EBk_shifted > 50 && sx3EUp > 50 && sx3EDn > 50) {
printf("sx3EUp: %f, sx3EDn: %f, sx3EBk_shifted: %f\n", sx3EUp, sx3EDn, sx3EBk_shifted);
printf("Filling hsx3uVsx3d_23_shifted: %f\n", sx3EUp / sx3EBk_shifted);
// hsx3uVsx3d_23->Fill(sx3EUp / sx3EBk_shifted, (-0.924754*sx3EUp+0.916671) / sx3EBk_shifted);
}
}
}
/* if ((sx3ChUp == 1 && sx3ChDn == 0)) {
if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
if (sx3EBk != 0 && sx3EBk > 50 && sx3EUp > 50 && sx3EDn > 50) {
printf("sx3EUp: %f, sx3EDn: %f, sx3EBk: %f\n", sx3EUp, sx3EDn, sx3EBk);
printf("Filling hsx3uVsx3d_01: %f\n", sx3EUp / sx3EBk);
hsx3uVsx3d_45->Fill(sx3EUp / sx3EBk, sx3EDn / sx3EBk);
}
}
}
else if ((sx3ChUp == 3 && sx3ChDn == 2)) {
if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
if (sx3EBk != 0 && sx3EBk > 50 && sx3EUp > 50 && sx3EDn > 50) {
printf("sx3EUp: %f, sx3EDn: %f, sx3EBk: %f\n", sx3EUp, sx3EDn, sx3EBk);
printf("Filling hsx3uVsx3d_23: %f\n", sx3EUp / sx3EBk);
hsx3uVsx3d_67->Fill(sx3EUp / sx3EBk, sx3EDn / sx3EBk);
}
}
}*/
if (sx3ChUp == 1 && sx3ChDn == 0){
//if (sx3ChUp == 1 || sx3ChDn == 0 || sx3ChUp == 3 || sx3ChDn == 2 || sx3ChUp == 5 || sx3ChDn == 4 || sx3ChUp == 7 || sx3ChDn == 6) {
if (sx3ChUp != -1 && sx3ChBk !=-1 && sx3ChDn !=-1) {
if (sx3EBk_shifted > 50 && sx3EUp > 50 && sx3EDn>50 &&sx3E_u_matched_01>50 && sx3E_u_matched_01>50) {
//printf("sx3EUp: %f, sx3EDn: %f, sx3E_u_matched_01: %f,sx3E_d_matched_01: %f\n", sx3EUp, sx3EDn, sx3E_u_matched_01,sx3E_d_matched_01);
printf("Filling hsx3uVsx3d_nn: %f, gggggg: %f \n", (sx3EUp+sx3EDn),(sx3E_u_matched_01+sx3E_d_matched_01) );
hsx3uVsx3d_45->Fill((sx3EUp+sx3E_d_matched_01),sx3EBk_shifted);
hsx3uVsx3d_67->Fill((sx3EUp+sx3E_d_matched_01),sx3E_fb_matched_01);
}
}
}
/*if (sx3ChBk > 8) {
if ((sx3ChUp == 7 && sx3ChDn == 6) ||
(sx3ChUp == 5 && sx3ChDn == 4) ||
(sx3ChUp == 3 && sx3ChDn == 2) ||
(sx3ChUp == 1 && sx3ChDn == 0)) {
if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
if (sx3EBk != 0 && sx3EBk > 50 && sx3EUp > 50 && sx3EDn > 50) {
hsx3uVsx3d->Fill(sx3EUp / sx3EBk, sx3EDn / sx3EBk);
hsx3Vsx3->Fill(sx3EUp ,sx3EDn);
printf("sx3EUp: %f | sx3EDn: %f | sx3EBk: %f | sx3ChUp: %d | sx3ChDn: %d | sx3ChBk: %d\n", sx3EUp, sx3EDn, sx3EBk, sx3ChUp, sx3ChDn, sx3ChBk);
}
}
}
}*/
//else {
//printf("sx3EUp\n");
//}
if (sx3ChUp == 1 && sx3ChDn == 0){
if (sx3ChUp != -1 && sx3ChBk !=-1 && sx3ChDn !=-1) {
if (sx3E_d_matched_01> sx3EUp ) {
//printf("hZd_01_1_dn: %f\n", sx3E_d_matched_01);
//printf("hZd_01_1_b: %f\n", sx3E_fb_matched_01);
hZd_01_1->Fill((2*(sx3E_d_matched_01+(coeff*diff))/sx3E_fb_matched_01)-1);
}
else if(sx3EUp> sx3E_d_matched_01) {
//printf("hZd_01_2_sx3EUp: %f\n",sx3EUp );
//printf("hZd_01_2_sx3EDn: %f\n",sx3E_fb_matched_01);
hZd_01_2->Fill(1-(2*(sx3EUp+(1-coeff)*diff))/sx3E_fb_matched_01);
}
else if(sx3EUp>0.0 && sx3E_d_matched_01>0.0 && sx3E_d_matched_01>=sx3EUp ) {
hZd_01_3->Fill((2*(sx3E_d_matched_01+ coeff*diff)/sx3E_fb_matched_01)-1);
}
else if(sx3EUp>0.0 && sx3E_d_matched_01>0.0 && sx3E_d_matched_01<sx3EUp ) {
hZd_01_4->Fill(1-(2*(sx3EUp+ (1-coeff)*diff)/sx3E_fb_matched_01));
}
}
}
}
for (int j = 0; j < pc.multi; j++) {
if (sx3.ch[index] > 8) {
hsx3VpcE->Fill(sx3.e[i], pc.e[j]);
}
}
}
sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
hitPos = sx3_contr.GetHitPos();
HitNonZero = true;
// hitPos.Print();
}
}
// //======================= PC
//########################################################### Track constrcution
//############################## DO THE KINEMATICS
}
return kTRUE;
}
void gainmatch::Terminate(){
}

114
gainmatch.h Normal file
View File

@ -0,0 +1,114 @@
#ifndef gainmatch_h
#define gainmatch_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include "Armory/ClassDet.h"
class gainmatch : public TSelector {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
// Fixed size dimensions of array or collections stored in the TTree if any.
// Declaration of leaf types
Det sx3;
Det qqq;
Det pc ;
ULong64_t evID;
UInt_t run;
// List of branches
TBranch *b_eventID; //!
TBranch *b_run; //!
TBranch *b_sx3Multi; //!
TBranch *b_sx3ID; //!
TBranch *b_sx3Ch; //!
TBranch *b_sx3E; //!
TBranch *b_sx3T; //!
TBranch *b_qqqMulti; //!
TBranch *b_qqqID; //!
TBranch *b_qqqCh; //!
TBranch *b_qqqE; //!
TBranch *b_qqqT; //!
TBranch *b_pcMulti; //!
TBranch *b_pcID; //!
TBranch *b_pcCh; //!
TBranch *b_pcE; //!
TBranch *b_pcT; //!
gainmatch(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~gainmatch() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree);
virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
virtual void SetOption(const char *option) { fOption = option; }
virtual void SetObject(TObject *obj) { fObject = obj; }
virtual void SetInputList(TList *input) { fInput = input; }
virtual TList *GetOutputList() const { return fOutput; }
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(gainmatch,0);
};
#endif
#ifdef gainmatch_cxx
void gainmatch::Init(TTree *tree){
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_eventID);
fChain->SetBranchAddress("run", &run, &b_run);
sx3.SetDetDimension(24,12);
qqq.SetDetDimension(4,32);
pc.SetDetDimension(2,24);
fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID);
fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch);
fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E);
fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T);
fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti);
fChain->SetBranchAddress("qqqID", &qqq.id, &b_qqqID);
fChain->SetBranchAddress("qqqCh", &qqq.ch, &b_qqqCh);
fChain->SetBranchAddress("qqqE", &qqq.e, &b_qqqE);
fChain->SetBranchAddress("qqqT", &qqq.t, &b_qqqT);
fChain->SetBranchAddress("pcMulti", &pc.multi, &b_pcMulti);
fChain->SetBranchAddress("pcID", &pc.id, &b_pcID);
fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh);
fChain->SetBranchAddress("pcE", &pc.e, &b_pcE);
fChain->SetBranchAddress("pcT", &pc.t, &b_pcT);
}
Bool_t gainmatch::Notify(){
return kTRUE;
}
void gainmatch::SlaveBegin(TTree * /*tree*/){
TString option = GetOption();
}
void gainmatch::SlaveTerminate(){
}
#endif // #ifdef gainmatch_cxx