diff --git a/Analysis.C b/Analysis.C index d646366..b9b575a 100644 --- a/Analysis.C +++ b/Analysis.C @@ -16,5 +16,5 @@ void Analysis(int start, int end) { // Define a macro with the same name as the script void Analysis() { - Analysis(72, 194); // Adjust the range if needed + Analysis(150, 194); // Adjust the range if needed } \ No newline at end of file diff --git a/Analyzer.C b/Analyzer.C index 7394293..df60a97 100644 --- a/Analyzer.C +++ b/Analyzer.C @@ -58,18 +58,17 @@ void Analyzer::Begin(TTree * /*tree*/){ hqqqVpcIndex->SetNdivisions( -612, "x"); 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( -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( -12, "y"); 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); - sx3_contr.ConstructGeo(); pw_contr.ConstructGeo(); @@ -283,7 +282,6 @@ Bool_t Analyzer::Process(Long64_t entry){ float aE = 0; float cE = 0; - bool multi_an =false; // if( ID[0].first < 1 ) { // aID = pc.ch[ID[0].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); - // for( int k = 0; k < qqq.multi; k++){ - // if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){ + for( int k = 0; k < qqq.multi; k++){ + if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){ + + int multi_an =0; for(int l=0;l=1){ + for(int l=0;l24){ + cE = E[l].second; + } + } + } + + } + } hanVScatsum->Fill(aE,cE); if( ID[0].first < 1 ) { diff --git a/Armory/#ClassPW.h# b/Armory/#ClassPW.h# new file mode 100644 index 0000000..ddfa498 --- /dev/null +++ b/Armory/#ClassPW.h# @@ -0,0 +1,283 @@ +#ifndef ClassPW_h +#define ClassPW_h + +#include +#include +#include + +struct PWHitInfo{ + std::pair nearestWire; // anode, cathode + std::pair nearestDist; // anode, cathode + + std::pair nextNearestWire; // anode, cathode + std::pair 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 GetNearestID() const {return hitInfo.nearestWire;} + std::pair GetNearestDistance() const {return hitInfo.nearestDist;} + std::pair Get2ndNearestID() const {return hitInfo.nextNearestWire;} + std::pair 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> An; // the anode wire position vector in space + std::vector> 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 p1; // anode + std::pair 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 \ No newline at end of file diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index 2cd6a38..de4fc58 100644 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -4,6 +4,7 @@ #include #include #include +#include struct PWHitInfo{ std::pair nearestWire; // anode, cathode @@ -65,8 +66,6 @@ public: 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); - double CircularMean(std::vector> wireList); - void Print(){ printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first, hitInfo.nearestDist.first, @@ -282,21 +281,4 @@ inline double PW::GetZ0(){ } -inline double PW::CircularMean(std::vector> 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 \ No newline at end of file diff --git a/TrackRecon.C b/TrackRecon.C index ebfe1b5..1f7ebae 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -1,6 +1,6 @@ -#define Analyzer_cxx +#define TrackRecon_cxx -#include "Analyzer.h" +#include "TrackRecon.h" #include #include #include @@ -14,20 +14,6 @@ #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; SX3 sx3_contr; @@ -37,7 +23,7 @@ bool HitNonZero; TH1F * hZProj; -void Analyzer::Begin(TTree * /*tree*/){ +void TrackRecon::Begin(TTree * /*tree*/){ TString option = GetOption(); hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600); @@ -47,17 +33,14 @@ void Analyzer::Begin(TTree * /*tree*/){ } - - - -Bool_t Analyzer::Process(Long64_t entry){ +Bool_t TrackRecon::Process(Long64_t entry){ // if ( entry > 100 ) return kTRUE; hitPos.Clear(); HitNonZero = false; - // if( entry > 1) return kTRUE; + if( entry > 1) return kTRUE; // printf("################### ev : %llu \n", entry); b_sx3Multi->GetEntry(entry); @@ -133,13 +116,6 @@ Bool_t Analyzer::Process(Long64_t entry){ }else{ 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); @@ -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 rho = 10.+40./16.*(chRing+0.5); // if(qqq.e[i]>50){ - hqqqPolar->Fill( theta, rho); + // hqqqPolar->Fill( theta, rho); // } // qqq.used[i] = true; // qqq.used[j] = true; @@ -187,7 +163,6 @@ Bool_t Analyzer::Process(Long64_t entry){ } // //======================= PC - PCHit_1An hitInfo; ID.clear(); int counter=0; @@ -197,33 +172,51 @@ Bool_t Analyzer::Process(Long64_t entry){ if( E.size()==3 ){ float aE = 0; float cE = 0; - bool multi_an =false; - for(int l=0;lFill(aE,cE); + if(multi_an==1){ + for(int l=0;l24){ + cE = E[l].second; + } + } + } + + //using CalTrack3 to get the track position and direction + + // hanVScatsum->Fill(aE,cE); if( HitNonZero){ - pw_contr.CalTrack3( hitPos, hitinfo, cID); - hZProj->Fill(pw_contr.GetZ0()); + if (ID.size() == 3) { + 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; } -void Analyzer::Terminate(){ +void TrackRecon::Terminate(){ gStyle->SetOptStat("neiou"); - TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000); - canvas->Divide(3,3); - - //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"); - + TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 200, 200); + padID=1; canvas->cd(padID); canvas->cd(padID)->SetGrid(1); + hZProj->Draw(); - // hanVScatsum->Draw("colz"); } diff --git a/TrackRecon.h b/TrackRecon.h new file mode 100644 index 0000000..3a93782 --- /dev/null +++ b/TrackRecon.h @@ -0,0 +1,114 @@ +#ifndef TrackRecon_h +#define TrackRecon_h + +#include +#include +#include +#include + +#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 diff --git a/archivist b/archivist new file mode 100644 index 0000000..4c800ef --- /dev/null +++ b/archivist @@ -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." diff --git a/gainmatch.C b/gainmatch.C new file mode 100644 index 0000000..8772d51 --- /dev/null +++ b/gainmatch.C @@ -0,0 +1,541 @@ +#define gainmatch_cxx + +#include "gainmatch.h" +#include +#include +#include +#include +#include + + +#include +#include + +#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> ID; // first = id, 2nd = index +for (int i = 0; i < sx3.multi; i++) { + ID.push_back(std::pair(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 &a, const std::pair &b) { + return a.first < b.first; + }); + + std::vector> 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_01Fill(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(){ + + +} \ No newline at end of file diff --git a/gainmatch.h b/gainmatch.h new file mode 100644 index 0000000..5d5697b --- /dev/null +++ b/gainmatch.h @@ -0,0 +1,114 @@ +#ifndef gainmatch_h +#define gainmatch_h + +#include +#include +#include +#include + +#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 \ No newline at end of file