ANASEN_analysis/TrackRecon.C

246 lines
5.8 KiB
C

#define TrackRecon_cxx
#include "TrackRecon.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <utility>
#include <algorithm>
#include "Armory/ClassSX3.h"
#include "Armory/ClassPC1An.h"
#include "TVector3.h"
int padID = 0;
SX3 sx3_contr;
PC pw_contr;
TVector3 hitPos;
bool HitNonZero;
TH1F * hZProj;
void TrackRecon::Begin(TTree * /*tree*/){
TString option = GetOption();
hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
sx3_contr.ConstructGeo();
pw_contr.ConstructGeo();
}
Bool_t TrackRecon::Process(Long64_t entry){
// if ( entry > 100 ) return kTRUE;
hitPos.Clear();
HitNonZero = 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
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;
} );
// printf("##############################\n");
// for( size_t i = 0; i < ID.size(); i++) printf("%zu | %d %d \n", i, ID[i].first, ID[i].second );
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]);
}
}
}
// printf("---------- sx3ID Multi : %zu \n", sx3ID.size());
if( found ){
int sx3ChUp, sx3ChDn, sx3ChBk;
float sx3EUp, sx3EDn;
// printf("------ sx3 ID : %d, multi: %zu\n", sx3ID[0].first, sx3ID.size());
for( size_t i = 0; i < sx3ID.size(); i++ ){
int index = sx3ID[i].second;
// printf(" %zu | index %d | ch : %d, energy : %d \n", i, index, sx3.ch[index], sx3.e[index]);
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];
}
}
sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
hitPos = sx3_contr.GetHitPos();
HitNonZero = true;
// hitPos.Print();
}
}
// //======================= QQQ
for( int i = 0; i < qqq.multi; i ++){
for( int j = i + 1; j < qqq.multi; j++){
if( qqq.id[i] == qqq.id[j] ){ // must be same detector
int chWedge = -1;
int chRing = -1;
if( qqq.ch[i] < qqq.ch[j]){
chRing = qqq.ch[j] - 16;
chWedge = qqq.ch[i];
}else{
chRing = qqq.ch[i];
chWedge = qqq.ch[j] - 16;
}
// printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing);
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);
// }
// qqq.used[i] = true;
// qqq.used[j] = true;
if( !HitNonZero ){
double x = rho * TMath::Cos(theta);
double y = rho * TMath::Sin(theta);
hitPos.SetXYZ(x, y, 23 + 75 + 30);
HitNonZero = true;
}
}
}
}
// //======================= PC
ID.clear();
int counter=0;
std::vector<std::pair<int, double>> E;
E.clear();
if( E.size()==3 ){
float aE = 0;
float cE = 0;
int multi_an =0;
for(int l=0;l<E.size();l++){
if(E[l].first<24 && E[l].first!=20 && E[l].first!=12){
multi_an++;
}
}
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 (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());
}
}
}
// }
}
//########################################################### Track constrcution
//############################## DO THE KINEMATICS
return kTRUE;
}
void TrackRecon::Terminate(){
gStyle->SetOptStat("neiou");
TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 200, 200);
padID=1;
canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
hZProj->Draw();
}