Commit 27020b00 authored by Jérémie Dudouet's avatar Jérémie Dudouet
Browse files

Update TreeBuilder to include PSA and CRYSTAL frames at local level in root trees

parent af3bf92d
......@@ -27,10 +27,201 @@
#include "TrackedFrame.h"
#include "PSAFrame.h"
#include "TGeoMatrix.h"
using namespace AGAPRO;
/// ****** CRYSTAL ******///
TB_AGATA_CRYSTAL::TB_AGATA_CRYSTAL(TString name, TString ADFKey) : TB_Detector(name,ADFKey)
{
gActualClass = "TB_AGATA_CRYSTAL";
}
TB_AGATA_CRYSTAL::~TB_AGATA_CRYSTAL()
{
}
void TB_AGATA_CRYSTAL::Init()
{
}
void TB_AGATA_CRYSTAL::InitValues()
{
Log.ClearMessage();
Log.GetProcessName() = gActualClass;
Log.SetProcessMethod("InitValues()");
crystalTS = 0;
crystalId = 0;
memset(SegmentTraces, 0, sizeof(SegmentTraces));
memset(CoreTraces, 0, sizeof(CoreTraces));
memset(SegE, 0, sizeof(SegE));
memset(CoreE, 0, sizeof(CoreE));
}
void TB_AGATA_CRYSTAL::InitTree()
{
Log.ClearMessage();
Log.GetProcessName() = gActualClass;
Log.SetProcessMethod("InitTree()");
if(fTree == nullptr)
return;
fTree->Branch("segTraces", SegmentTraces, Form("segTraces[%d]/F",kNSG*kNSamples));
fTree->Branch("coreTraces", CoreTraces, Form("coreTraces[%d]/F",kNCC*kNSamples));
fTree->Branch("segE", SegE, Form("segE[%d]/F",kNSG));
fTree->Branch("coreE", CoreE, Form("coreE[%d]/F",kNCC));
fTree->Branch("crystalId", &crystalId, "crystalId/I");
fTree->Branch("crystalTS", &crystalTS, "crystalTS/l");
}
void TB_AGATA_CRYSTAL::Process(Int_t idet)
{
Log.ClearMessage();
Log.GetProcessName() = "TB_AGATA_CRYSTAL";
Log.SetProcessMethod(Form("Process(%d)",idet));
fSFrame = fTrigger->GetInputSharedFP();
fSFrame->GetFrame()->Read();
CrystalInterface *data = GetDataPointer<CrystalInterface>(fSFrame);
crystalTS = ((AgataKey *)fSFrame->GetFrame()->GetKey())->GetTimeStamp();
crystalId = data->GetUID();
if(!WarnedDone && data->GetCore(0)->GetSignal()->GetLength()!=kNSamples) {
Log << error << " Signal length: " << data->GetCore(0)->GetSignal()->GetLength() << " ; Stored length in tree: " << kNSamples << nline;
WarnedDone = true;
}
for (Int_t i = 0; i < kNCC; i++ ) {
memcpy(&CoreTraces[i*kNSamples],data->GetCore(i)->GetSignal()->Address(CoreTraces),sizeof(Float_t)*kNSamples);
CoreE[i] = data->GetCore(i)->GetE();
}
for (Int_t i = 0; i < kNSG; i++ ) {
memcpy(&SegmentTraces[i*kNSamples],data->GetSegment(i)->GetSignal()->Address(SegmentTraces),sizeof(Float_t)*kNSamples);
SegE[i] = data->GetSegment(i)->GetE();
}
}
/// ****** PSA ******///
TB_AGATA_PSA::TB_AGATA_PSA(TString name, TString ADFKey) : TB_Detector(name,ADFKey)
{
gActualClass = "TB_AGATA_PSA";
}
TB_AGATA_PSA::~TB_AGATA_PSA()
{
}
void TB_AGATA_PSA::Init()
{
ReadTransformation();
}
void TB_AGATA_PSA::InitValues()
{
Log.ClearMessage();
Log.GetProcessName() = gActualClass;
Log.SetProcessMethod("InitValues()");
nbHits = 0;
coreId = 0.;
coreE0 = 0.;
coreE1 = 0.;
coreT0 = 0.;
coreT1 = 0.;
coreTS = 0.;
}
void TB_AGATA_PSA::InitTree()
{
Log.ClearMessage();
Log.GetProcessName() = gActualClass;
Log.SetProcessMethod("InitPSAHits()");
if(fTree == nullptr)
return;
fTree->Branch("nbHits", &nbHits, "nbHits/I");
fTree->Branch("hitE", hitE, "hitE[nbHits]/F");
fTree->Branch("hitX", hitX, "hitX[nbHits]/F");
fTree->Branch("hitY", hitY, "hitY[nbHits]/F");
fTree->Branch("hitZ", hitZ, "hitZ[nbHits]/F");
fTree->Branch("hitGX", hitGX, "hitGX[nbHits]/F");
fTree->Branch("hitGY", hitGY, "hitGY[nbHits]/F");
fTree->Branch("hitGZ", hitGZ, "hitGZ[nbHits]/F");
fTree->Branch("hitSg", hitSg, "hitSg[nbHits]/I");
fTree->Branch("coreId", &coreId, "coreId/I");
fTree->Branch("coreE0", &coreE0, "coreE0/F");
fTree->Branch("coreE1", &coreE1, "coreE1/F");
fTree->Branch("coreT0", &coreT0, "coreT0/F");
fTree->Branch("coreT1", &coreT1, "coreT1/F");
fTree->Branch("coreTS", &coreTS, "coreTS/l");
}
void TB_AGATA_PSA::Process(Int_t idet)
{
Log.ClearMessage();
Log.GetProcessName() = "TB_AGATA_PSA";
Log.SetProcessMethod(Form("Process(%d)",idet));
fSFrame = fTrigger->GetInputSharedFP();
fSFrame->GetFrame()->Read();
PSAInterface *data = GetDataPointer<PSAInterface>(fSFrame);
nbHits = data->GetNbHits();
if ( nbHits > MaxHits ) {
Log << error << " Number of Hits : " << nbHits << " ; Max nbHits = " << MaxHits << nline;
InitValues();
return;
}
// Core like branches
coreId = data->GetUID();
coreE0 = data->GetE(0);
coreE1 = data->GetE(1);
coreT0 = data->GetT(0);
coreT1 = data->GetT(1);
coreTS = ((AgataKey *)fSFrame->GetFrame()->GetKey())->GetTimeStamp();
coreTS = ((AgataKey *)fSFrame->GetFrame()->GetKey())->GetTimeStamp();
// Loop on the hits
for (Int_t j = 0; j < nbHits ; ++j) {
const PSAHit* pHit = (PSAHit*)data->GetHit(j);
hitE[j] = pHit->GetE();
hitX[j] = pHit->GetX();
hitY[j] = pHit->GetY();
hitZ[j] = pHit->GetZ();
hitSg[j] = pHit->GetID();
//Get hit position
Double_t xLocal, yLocal, zLocal, xGlobal, yGlobal, zGlobal;
pHit->GetXYZ(xLocal, yLocal, zLocal);
//To the AGATA global ref
Local2Global(data->GetUID(), xLocal, yLocal, zLocal, xGlobal, yGlobal, zGlobal);
// Set the Position of the hit in global ref
hitGX[j]= xGlobal; // inversion X-Y and opposite direction
hitGY[j]= yGlobal; // inversion Y-X and opposite direction
hitGZ[j]= zGlobal; // z
}
}
/// ****** Builder ******///
TB_AGATA_Builder::TB_AGATA_Builder(TString name, TString ADFKey) : TB_Detector(name,ADFKey)
......@@ -40,13 +231,11 @@ TB_AGATA_Builder::TB_AGATA_Builder(TString name, TString ADFKey) : TB_Detector(n
TB_AGATA_Builder::~TB_AGATA_Builder()
{
if(fMatrixList)
delete fMatrixList;
}
void TB_AGATA_Builder::Init()
{
fFrame = fTrigger->AddUtility("data:psa",ConfAgent::theGlobalAgent());
fSFrame = fTrigger->AddUtility("data:psa",ConfAgent::theGlobalAgent());
ReadTransformation();
}
......@@ -150,12 +339,9 @@ void TB_AGATA_Builder::InitTree()
temparray[43] = "41 44 42";
temparray[44] = "13 16 42 43 41 14";
for(int i=0 ; i<45 ; i++)
{
for(int i=0 ; i<45 ; i++) {
TObjArray *arr = temparray[i].ReplaceAll("\t"," ").Tokenize(" ");
for(int j=0 ; j<arr->GetEntries() ; j++)
{
for(int j=0 ; j<arr->GetEntries() ; j++) {
Int_t det = ((TString)arr->At(j)->GetName()).Atoi();
TString name = Form("%.2d:%.2d",i,det);
fAreNeighbours[name] = true;
......@@ -168,9 +354,9 @@ void TB_AGATA_Builder::Process(Int_t idet)
{
Log.ClearMessage();
Log.GetProcessName() = "TB_AGATA_Builder";
Log.SetProcessMethod("Process()");
Log.SetProcessMethod(Form("Process(%d)",idet));
if(fNoMergerMode)
if(!fMergerMode)
fCompFrame = (AgataCompositeFrame*) fTrigger->GetInputSharedFP()->GetFrame();
else
fCompFrame = (AgataCompositeFrame*) fTrigger->GetInputSharedFP(idet+1)->GetFrame();
......@@ -183,8 +369,6 @@ void TB_AGATA_Builder::Process(Int_t idet)
nbCores = fCompFrame->GetNbSubFrame();
// cout<<"builder "<<TSHit<<" "<<nbCores<<endl;
Float_t MaxE[nbCores];
Float_t coreX[nbCores];
Float_t coreY[nbCores];
......@@ -195,17 +379,15 @@ void TB_AGATA_Builder::Process(Int_t idet)
memset(coreY, 0., sizeof(coreY));
memset(coreZ, 0., sizeof(coreZ));
for(Int_t i=0 ; i < nbCores ; i++)
{
fCompFrame->LinkSubFrame(i,fFrame->GetFrame());
fFrame->GetFrame()->Read();
for(Int_t i=0 ; i < nbCores ; i++) {
fCompFrame->LinkSubFrame(i,fSFrame->GetFrame());
fSFrame->GetFrame()->Read();
PSAInterface *data = GetDataPointer<PSAInterface>(fFrame);
PSAInterface *data = GetDataPointer<PSAInterface>(fSFrame);
nbHitsperCry[i] = data->GetNbHits();
if ( nbHitsperCry[i] == MaxHits )
{
if ( nbHitsperCry[i] > MaxHits ) {
Log << error << " Number of Hits to hit in : " << nbHitsperCry[i] << "Max nbHits = " << MaxHits << nline;
InitValues();
......@@ -217,11 +399,10 @@ void TB_AGATA_Builder::Process(Int_t idet)
coreE0[i] = data->GetE(0);
coreE1[i] = data->GetE(1);
coreT0[i] = data->GetT(0);
coreTS[i] = ((AgataKey *)fFrame->GetFrame()->GetKey())->GetTimeStamp();
coreTS[i] = ((AgataKey *)fSFrame->GetFrame()->GetKey())->GetTimeStamp();
// Loop on the hit to get over the hits with the max of the energy
for (Int_t j = 0; j < nbHitsperCry[i]; ++j)
{
for (Int_t j = 0; j < nbHitsperCry[i]; ++j) {
const PSAHit* pHit = (PSAHit*)data->GetHit(j);
hitE[nbHits+j] = pHit->GetE();
......@@ -236,15 +417,14 @@ void TB_AGATA_Builder::Process(Int_t idet)
//To the AGATA global ref
Local2Global(data->GetUID(), xLocal, yLocal, zLocal, xGlobal, yGlobal, zGlobal);
if(pHit->GetE()>MaxE[i])
{
MaxE[i] = pHit->GetE();
coreX[i] = xGlobal;
coreY[i] = yGlobal;
coreZ[i] = zGlobal;
if(pHit->GetE()>MaxE[i]) {
MaxE[i] = pHit->GetE();
coreX[i] = xGlobal;
coreY[i] = yGlobal;
coreZ[i] = zGlobal;
}
// Set the Position of the hit in global ref of PRESPEC
// Set the Position of the hit in global ref
hitGX[nbHits+j]= xGlobal; // inversion X-Y and opposite direction
hitGY[nbHits+j]= yGlobal; // inversion Y-X and opposite direction
hitGZ[nbHits+j]= zGlobal; // z
......@@ -254,10 +434,8 @@ void TB_AGATA_Builder::Process(Int_t idet)
}
//Addback calcluation
std::vector< core_property > corelist;
for(int i=0 ; i<nbCores ; i++)
{
for(int i=0 ; i<nbCores ; i++) {
core_property acore;
acore.E = coreE0[i];
acore.Id = coreId[i];
......@@ -273,16 +451,12 @@ void TB_AGATA_Builder::Process(Int_t idet)
std::sort(corelist.begin(), corelist.end(), coreCompare);
for(int i=0 ; i<nbCores ; i++)
{
for(int j=i+1 ; j<nbCores ; j++)
{
if(fAreNeighbours[Form("%.2d:%.2d",corelist[i].Id,corelist[j].Id)])
{
for(int i=0 ; i<nbCores ; i++) {
for(int j=i+1 ; j<nbCores ; j++) {
if(fAreNeighbours[Form("%.2d:%.2d",corelist[i].Id,corelist[j].Id)]) {
corelist[i].E += corelist[j].E;
if(corelist[j].HighHitE > corelist[i].HighHitE)
{
if(corelist[j].HighHitE > corelist[i].HighHitE) {
corelist[i].HighHitE = corelist[j].HighHitE;
corelist[i].HighHitX = corelist[j].HighHitX;
corelist[i].HighHitY = corelist[j].HighHitY;
......@@ -298,10 +472,8 @@ void TB_AGATA_Builder::Process(Int_t idet)
std::sort(corelist.begin(), corelist.end(), coreCompare);
}
for(int i=0 ; i<nbCores ; i++)
{
if(corelist[i].E>0)
{
for(int i=0 ; i<nbCores ; i++) {
if(corelist[i].E>0) {
AddE[nbAdd] = corelist[i].E;
AddX[nbAdd] = corelist[i].HighHitX;
AddY[nbAdd] = corelist[i].HighHitY;
......@@ -314,126 +486,6 @@ void TB_AGATA_Builder::Process(Int_t idet)
}
}
TGeoCombiTrans* TB_AGATA_Builder::FillTransMatrix(Double_t* trans, Double_t* rot)
{
TGeoRotation rotation; rotation.SetMatrix(rot);
TGeoTranslation translation;
translation.SetTranslation(trans[0], trans[1], trans[2]);
return new TGeoCombiTrans(translation, rotation);
}
void TB_AGATA_Builder::Local2Global(Int_t detID,
Double_t xl, Double_t yl, Double_t zl,
Double_t& xg, Double_t& yg, Double_t& zg)
{
Log.ClearMessage();
Log.GetProcessName() = gActualClass;
Log.SetProcessMethod("Local2Global()");
if (detID < 0 || detID > 180)
{
Log << error << "Wrong detector id number" << dolog;
return;
}
TGeoCombiTrans* mat = static_cast<TGeoCombiTrans*> ( fMatrixList->At(detID) );
Double_t local[3] = {xl, yl, zl};
Double_t global[3] = {0., 0., 0.};
mat->LocalToMaster(local, global);
xg = global[0];
yg = global[1];
zg = global[2];
}
bool TB_AGATA_Builder::ReadTransformation()
{
Log.ClearMessage();
Log.GetProcessName() = gActualClass;
Log.SetProcessMethod("ReadTransformation()");
Log << info << " Ge positions extracted from CrystalPositionLookUpTable " << nline;
TString FileName = fConfPath + "/CrystalPositionLookUpTable";
ifstream in(FileName, ios::in);
if (!in) {
Log << error << Form("file %s not found\n", FileName.Data()) << dolog;
return false;
}
fMatrixList = new TObjArray(180);
fMatrixList->SetOwner();
Int_t nd, k;
Double_t x,y,z;
Double_t trans[3];
Double_t rot[9];
Char_t line[255];
while ( !in.eof() ) {
in.getline(line,255);
if((((TString)line)==""))
continue;
if(sscanf(line, "%d %d %lf %lf %lf", &nd, &k, &x, &y, &z) != 5) {
Log << error << "a Problem reading translation factors" << dolog;
cin.get();
return false;
}
trans[0] = x; trans[1] = y, trans[2] = z;
if (nd > 180) {
Log << error << Form("Id %d too large, max: %d\n", nd, 180) << dolog;
return false;
}
in.getline(line,255);
if(sscanf(line, "%d %lf %lf %lf", &k, &x, &y, &z) != 4)
{
Log << error << Form("b Problem reading 1st rotation factors for id %d\n", nd) << dolog;
return false;
}
rot[0] = x; rot[1] = y; rot[2] = z;
in.getline(line,255);
if(sscanf(line, "%d %lf %lf %lf", &k, &x, &y, &z) != 4)
{
Log << error << Form("c Problem reading 2nd rotation factors for id %d\n", nd) << dolog;
return false;
}
rot[3] = x; rot[4] = y; rot[5] = z;
in.getline(line,255);
if(sscanf(line, "%d %lf %lf %lf", &k, &x, &y, &z) != 4)
{
Log << error << Form("d Problem reading 3rd rotation factors for id %d\n", nd) << dolog;
return false;
}
rot[6] = x; rot[7] = y; rot[8] = z;
// add matrix for the given id
if (fMatrixList->At(nd)) {
Log << error << Form("Matrix transformation already exist for id %d\n", nd) << dolog;
return false;
} else
fMatrixList->AddAt( FillTransMatrix(trans, rot), nd );
}
in.close( );
Log << dolog;
}
/// ****** Tracking ******///
......@@ -449,7 +501,7 @@ TB_AGATA_Tracking::~TB_AGATA_Tracking()
void TB_AGATA_Tracking::Init()
{
fFrame = fTrigger->AddUtility("data:tracked",ConfAgent::theGlobalAgent());
fSFrame = fTrigger->AddUtility("data:tracked",ConfAgent::theGlobalAgent());
}
void TB_AGATA_Tracking::InitValues()
......@@ -492,39 +544,33 @@ void TB_AGATA_Tracking::Process(Int_t idet)
{
Log.ClearMessage();
Log.GetProcessName() = gActualClass;
Log.SetProcessMethod("Process()");
Log.SetProcessMethod(Form("Process(%d)",idet));
if(fNoMergerMode)
{
fFrame = fTrigger->GetInputSharedFP();
}
if(!fMergerMode)
fSFrame = fTrigger->GetInputSharedFP();
else
fFrame = fTrigger->GetInputSharedFP(idet+1);
fSFrame = fTrigger->GetInputSharedFP(idet+1);
fFrame->GetFrame()->Read();
fSFrame->GetFrame()->Read();
// to get the data part of the frame
const GammaTrackedInterface *data = GetCstDataPointer<GammaTrackedInterface>(fFrame);
const GammaTrackedInterface *data = GetCstDataPointer<GammaTrackedInterface>(fSFrame);
if(data == nullptr)
return;
nbTrack = data->GetNbGamma();
if(nbTrack>MaxTrackedGamma)
{
if(nbTrack>MaxTrackedGamma) {
Log << error << " Number of tracked gamma : " << nbTrack << "Max nbTrack = " << MaxTrackedGamma << " event ignored !" << nline;
InitValues();
return;
}
//! get time stamp
TSTrack = ((AgataKey *)fFrame->GetFrame()->GetKey())->GetTimeStamp();
// cout<<"tracking "<< TSTrack << " " << nbTrack << endl;
TSTrack = ((AgataKey *)fSFrame->GetFrame()->GetKey())->GetTimeStamp();
for (UShort_t i = 0u; i < nbTrack; i++)
{
for (UShort_t i = 0u; i < nbTrack; i++) {
const TrackedHit *gamma1 = data->GetGamma(i);
const Hit *hit1 = gamma1->GetHit();
......@@ -546,7 +592,7 @@ void TB_AGATA_Tracking::Process(Int_t idet)
trackY2[i] = gamma1->GetY();
trackZ2[i] = gamma1->GetZ();
}
else{
else {
trackX2[i] = trackX1[i];
trackY2[i] = trackY1[i];
trackZ2[i] = trackZ1[i];
......@@ -557,8 +603,7 @@ void TB_AGATA_Tracking::Process(Int_t idet)
nbTrack--;
}
if(nbTrack==0)
{
if(nbTrack==0) {
InitValues();
return;
}
......
......@@ -21,12 +21,80 @@
#ifndef _TB_AGATA
#define _TB_AGATA
#include "AGAPRO_TB_Detector.h"
#include "CrystalFrame.h"
class TGeoCombiTrans;
#include "AGAPRO_TB_Detector.h"
namespace AGAPRO {
/// ****** CRYSTAL ******///
class TB_AGATA_CRYSTAL : public TB_Detector
{
protected:
static const int kNCC = ADF::CrystalInterface::kNbCores;
static const int kNSG = ADF::CrystalInterface::kNbSegments;
static const int kNSamples = 60;
private:
Float_t SegmentTraces[kNSG*kNSamples];
Float_t CoreTraces[kNCC*kNSamples];
Float_t SegE[kNSG];
Float_t CoreE[kNCC];
ULong64_t crystalId;
ULong64_t crystalTS;
bool WarnedDone = false;
public:
TB_AGATA_CRYSTAL(TString name="", TString ADFKey="");
virtual ~TB_AGATA_CRYSTAL();
void Init();
void InitValues();