#include "BrCombineTrack.h"
#include "BrDetectorVolume.h"
#include "BrEventNode.h"
#include "BrDataTable.h"
#include "BrDetectorTrack.h"
#include "BrTrack.h"
#include "BrTrackDefines.h"
ClassImp(BrCombineTrack)
///////////////////////////////////////////////////////////////////////
// //
// The BrCombineTrack class provides an interface to the BRAHMS //
// Tracking Combining between detectors //
// //
///////////////////////////////////////////////////////////////////////
BrCombineTrack::BrCombineTrack():BrModule()
{
// Default Constructor.
for(Int_t i=0;i<max_num_volumes;i++) fVolume_p[i] = NULL;
}
BrCombineTrack::BrCombineTrack(Text_t *Name,Char_t *Title)
:BrModule(Name,Title)
{
//Set volume parameters to null,
//Create clones arrays for internal use
//Read simulation parameters from default file (simpar.dat) in current directory
Int_t i;
for(i=0;i<max_num_volumes;i++) fVolume_p[i] = NULL;
CreateClonesArrays();
ReadSimulationParameters("simpar.dat"); //default filename to read
fFieldParametersRead = kFALSE;
}
BrCombineTrack::BrCombineTrack(Text_t *Name,Char_t *Title,Char_t *simparfile,Char_t *fieldfile)
:BrModule(Name,Title)
{
//Set volume parameters to null,
//Create clones arrays for internal use
//Read simulation parameters from file specified in parameter
Int_t i;
for(i=0;i<max_num_volumes;i++) fVolume_p[i] = NULL;
CreateClonesArrays();
ReadSimulationParameters(simparfile); //user specified filename to read
ReadFieldParameters(fieldfile); //user specified filename to read
}
BrCombineTrack::BrCombineTrack(Char_t *Name,Char_t *Title,Char_t *simparfile,Char_t *fieldfile,
BrDetectorVolume *tpc1,BrDetectorVolume *tpc2,
BrDetectorVolume *dc3,BrDetectorVolume *dc4,BrDetectorVolume *dc5,
BrDetectorVolume *d1,BrDetectorVolume *d2,BrDetectorVolume *d3,BrDetectorVolume *d4,
BrDetectorVolume *m0) : BrModule(Name,Title) {
//Set volume parameters to null,
//Create clones arrays for internal use
//Read simulation parameters from file specified in parameter
//Set volume parameters as given in arguments. This contstructor eliminates the need
//to use 10 statements of AddDetectorVolume
Int_t i;
for(i=0;i<max_num_volumes;i++) fVolume_p[i] = NULL;
CreateClonesArrays();
ReadSimulationParameters(simparfile); //user specified filename to read
ReadFieldParameters(fieldfile); //user specified filename to read
AddDetectorVolume(tpc1);
AddDetectorVolume(tpc2);
AddDetectorVolume(dc3);
AddDetectorVolume(dc4);
AddDetectorVolume(dc5);
AddDetectorVolume(d1);
AddDetectorVolume(d2);
AddDetectorVolume(d3);
AddDetectorVolume(d4);
AddDetectorVolume(m0);
}
void BrCombineTrack::CreateClonesArrays()
{
//Create the clones arrays used for internal structures.
Tracks12 = new BrClonesArray("BrTrack12",1000);
T2Match = new BrClonesArray("BrT2Match",1000);
Tracks34 = new BrClonesArray("BrTrack34",1000);
Tracks45 = new BrClonesArray("BrTrack45",1000);
Tracks345 = new BrClonesArray("BrTrack345",1000);
}
void BrCombineTrack::ClearClonesArrays()
{
//Clear the clones arrays used for internal structures.
Tracks12->Clear();
T2Match->Clear();
Tracks34->Clear();
Tracks45->Clear();
Tracks345->Clear();
}
BrCombineTrack::~BrCombineTrack()
{
//Destructor
//Delete the fVolume_p parameters
Int_t i;
if( fVolume_p != 0) {
for(i=0;i<max_num_volumes;i++) delete fVolume_p[i];
}
if(Tracks12) delete Tracks12;
if(Tracks34) delete Tracks34;
if(Tracks45) delete Tracks45;
if(T2Match) delete T2Match;
}
void BrCombineTrack::Event(BrEventNode* InputTable,BrEventNode* OutputTable)
{
//Execute the actions necessary for combining the different tracks.
Int_t i,volnotsetup,fspec;
static const Char_t *testvolname[] = {"T1","T2","T3","T4","T5","D1","D2","D3","D4","M0"};
BrDetectorVolume *testvol;
fspec = 1; //for a test for now...
ClearClonesArrays();
//First check whether volumenames are setup correctly (for that matter, set up at all!!)
volnotsetup = 0;
for(i=0;i<10;i++) {
testvol = GetDetectorVolume(testvolname[i]);
if(!testvol) {
cout<<"Detector Volume "<<testvolname[i]<<" is not setup yet"<<endl;
cout<<"Please set it up using AddDetectorVolume(...)"<<endl;
volnotsetup++;
}
}
//if(volnotsetup > 0) return;
if( fspec ) {
if( GetTrackingLevel("T12") == 1 ) {
CombineT1T2(InputTable);
Track12ToTarget();
}
if( GetTrackingLevel("T345") == 1 ) {
CombineT3T4T5(InputTable,OutputTable);
if( GetTrackingLevel("T2M") >= 1 ) {
CombineD12D34(InputTable,OutputTable);
Track345ToTarget();
}
}
}
else {
if( GetTrackingLevel("T12") == 1 ) {
CombineMIDST1T2(InputTable);
Track12MIDSToTarget();
}
}
}
void BrCombineTrack::CombineT1T2(BrEventNode *InputTable)
{
//Loop over the tracks found in T2 and convert the local to global coordinates,
//project the track to the center of the magnet.
//It will project the tracks in T1 forward to the center plane of
//the magnet. It will do the calculation in the local coordinate
//system of D2 which is most useful for this purpose.
//It will make simple corrections due to path length differences.
BrDetectorTrack* t1tr_p;
BrDetectorTrack* t2tr_p;
BrTrack12* track12_p;
Int_t i,j;
Float_t d2x,d2y,d2l,t;
Float_t xgt1[3],agt1[3];
Float_t xgt2[3],agt2[3];
#define max_t1 100
#define max_t2 50
Float_t xlt1[max_t1][3],alt1[max_t1][3];
Float_t xlt2[max_t2][3],alt2[max_t2][3];
Float_t xpt1[max_t1],ypt1[max_t1];
Float_t xpt2[max_t2],ypt2[max_t2];
Float_t cs1,cs2,sin1,sin2,dx,dy,aly,p;
Int_t status_t1t2;
Int_t it1,it2;
BrDetectorVolume *t1;
BrDetectorVolume *t2;
BrDetectorVolume *d2;
BrDataTable *T1DetectorTracks;
BrDataTable *T2DetectorTracks;
Int_t NumT1Tr,NumT2Tr;
T1DetectorTracks = InputTable->GetDataTable("T1 Detector Tracks");
if(T1DetectorTracks == 0) {
cout<<"Can't process CombineT1T2 because T1 Detector Tracks table not supplied"<<endl;
}
T2DetectorTracks = InputTable->GetDataTable("T2 Detector Tracks");
if(T2DetectorTracks == 0) {
cout<<"Can't process CombineT1T2 because T2 Detector Tracks table not supplied"<<endl;
}
if(!T1DetectorTracks||!T2DetectorTracks) return;
t1 = GetDetectorVolume("T1");
if(!t1) cout<<"Can't process CombineT1T2 because T1 volume not supplied"<<endl;
t2 = GetDetectorVolume("T2");
if(!t2) cout<<"Can't process CombineT1T2 because T2 volume not supplied"<<endl;
d2 = GetDetectorVolume("D2");
if(!d2) cout<<"Can't process CombineT1T2 because D2 volume not supplied"<<endl;
//Find the D2 sizes in the first time
d2x = d2->GetSize()[0] / (Float_t)2.0;
d2y = d2->GetSize()[1] / (Float_t)2.0;
d2l = d2->GetSize()[2] / (Float_t)2.0;
//Loop over tracks from t2
NumT1Tr = T1DetectorTracks->GetEntries();
NumT2Tr = T2DetectorTracks->GetEntries();
for(i=0;i<NumT2Tr;i++) {
t2tr_p = (BrDetectorTrack*)T2DetectorTracks->At(i);
t2->LocalToGlobal(t2tr_p->GetPos(),xgt2,0);
t2->LocalToGlobal(t2tr_p->GetAlph(),agt2,1);
// Convert into system for magnet D2
d2->GlobalToLocal(xgt2,xlt2[i],0);
d2->GlobalToLocal(agt2,alt2[i],1);
// Track to center of magnet D2
t = xlt2[i][2] / alt2[i][2];
xpt2[i] = xlt2[i][0] - alt2[i][0] * t;
ypt2[i] = xlt2[i][1] - alt2[i][1] * t;
}
//Loop over tracks from t1
for(i=0;i<NumT1Tr;i++) {
t1tr_p = (BrDetectorTrack*)T1DetectorTracks->At(i);
t1->LocalToGlobal(t1tr_p->GetPos(),xgt1,0);
t1->LocalToGlobal(t1tr_p->GetAlph(),agt1,1);
// Convert into system for magnet D2
d2->GlobalToLocal(xgt1,xlt1[i],0);
d2->GlobalToLocal(agt1,alt1[i],1);
t = -xlt1[i][2] / alt1[i][2];
xpt1[i] = xlt1[i][0] + alt1[i][0] * t;
ypt1[i] = xlt1[i][1] + alt1[i][1] * t;
}
//Loop over combinations of T1 and T2 tracks and look at their
//matching at the center. If the tracks are outside the volume
//they are rejected. Otherwise they are given a status/flag bit
for(it1=0;it1<NumT1Tr;it1++) {
t1tr_p = (BrDetectorTrack*)T1DetectorTracks->At(it1);
if(TMath::Abs(xpt1[it1])<d2x&&TMath::Abs(ypt1[it1])<d2y) {
// Loop over T2 tracks
for(it2=0;it2<NumT2Tr;it2++) {
t2tr_p = (BrDetectorTrack*)T2DetectorTracks->At(it2);
if(TMath::Abs(xpt2[it2])<d2x&&TMath::Abs(ypt2[it2])<d2y) {
dx = xpt2[it2] - xpt1[it1];
dy = ypt2[it2] - ypt1[it1];
aly = alt2[it2][1] - alt1[it1][1];
// Match tracks and calculate momentum
if(TMath::Abs(dx) < t1t2_par[0] &&
TMath::Abs(dy) < t1t2_par[1] &&
TMath::Abs(aly) < t1t2_par[3] ) {
// Set some status flag
status_t1t2 = 0;
if((TMath::Abs(xpt2[it2]) > (d2x-t1t2_par[0]))||
(TMath::Abs(ypt2[it2]) > (d2y-t1t2_par[1]))) {
status_t1t2 |= 1<<sf_t1_d2_midcut;
}
if((TMath::Abs(xpt1[it1]) > (d2x-t1t2_par[0]))||
(TMath::Abs(ypt1[it1]) > (d2y-t1t2_par[1]))) {
status_t1t2 |= 1<<sf_t2_d2_midcut;
}
cs1 = (Float_t)0.0;
cs2 = (Float_t)0.0;
for(j=0;j<3;j+=2) {
cs1 += alt1[it1][j] * alt1[it1][j];
cs2 += alt2[it2][j] * alt2[it2][j];
}
cs1 = (Float_t)TMath::Sqrt(cs1);
cs2 = (Float_t)TMath::Sqrt(cs2);
sin1 = -alt1[it1][0] / cs1;
sin2 = -alt2[it2][0] / cs2;
p = b2 * (d2l+d2l) / ((Float_t)3335.6*(sin1-sin2));
p = p / (Float_t)TMath::Sqrt(1.0-alt1[it1][1]*alt1[it1][1]);
BrClonesArray &track12 = *Tracks12;
Int_t ntrack12 = Tracks12->GetEntries();
track12_p = new(track12[ntrack12]) BrTrack12();
track12_p->SetP(p);
track12_p->SetDx(dx);
track12_p->SetDy(dy);
track12_p->SetT1tr(t1tr_p);
track12_p->SetT2tr(t2tr_p);
}
}
}
}
else {
// t1 track rejected, put appropriate stuff here
}
}
}
void BrCombineTrack::Track12ToTarget()
{
//This routine will loop over all tracks found using their position
//in T1 as a starting point for backward tracking. It will
//calculate Tx and Ty at the target position, calculate the
//scattering angle using the information for the track positions
//at the entrance to D1 ie on a virtual plane.
BrDetectorTrack* t1tr_p;
BrTrack12* track12_p;
Int_t NumTrack12;
Int_t it;
Float_t xl[3],al[3],xg[3],ag[3],aln[3];
Float_t d1x,d1y,d1l;
Float_t t,xc,yc,xe,ye,xt,yt,cs,si,cs2,sin1,sin2,r,th1,th2,theta,phi;
BrDetectorVolume *d1;
BrDetectorVolume *t1;
d1 = GetDetectorVolume("D1");
if(!d1) cout<<"Can't process TrackT12ToTarget because D1 volume not supplied"<<endl;
t1 = GetDetectorVolume("T1");
if(!t1) cout<<"Can't process TrackT12ToTarget because T1 volume not supplied"<<endl;
d1x = d1->GetSize()[0] / (Float_t)2.0;
d1y = d1->GetSize()[1] / (Float_t)2.0;
d1l = d1->GetSize()[2] / (Float_t)2.0;
//Loop over tracks
NumTrack12 = Tracks12->GetEntries();
for(it=0;it<NumTrack12;it++) {
track12_p = (BrTrack12*)Tracks12->At(it);
t1tr_p = track12_p->GetT1tr();
// The local track coordinates are given in the local coordinates
// for T1. Convert to global --> local for D1
t1->LocalToGlobal(t1tr_p->GetPos(),xg,0);
t1->LocalToGlobal(t1tr_p->GetAlph(),ag,1);
d1->GlobalToLocal(xg,xl,0);
d1->GlobalToLocal(ag,al,1);
// Project to front of d1 and bend the particle according to
// p and B.dl
t = - (xl[2] - d1l) / al[2];
xc = xl[0] + t * al[0];
yc = xl[1] + t * al[1];
cs2 = (Float_t)sqrt(al[0]*al[0] + al[2]*al[2]);
sin2 = al[0] / cs2;
r = (Float_t)3335.6 * track12_p->GetP() * cs2 / b1;
sin1 = -(Float_t)2.0 * d1l / r + sin2;
th1 = (Float_t)TMath::ASin(sin1);
th2 = (Float_t)TMath::ASin(sin2);
theta = th2 - th1;
// Rotate the coordinates for (al[0],al[2]) by the scattering angle
cs = (Float_t)TMath::Cos(theta);
si = (Float_t)TMath::Sin(theta);
aln[0] = al[0] * cs - al[2] * si;
aln[2] = al[1] * si - al[2] * cs;
aln[1] = al[1];
// These are now the direction cosines
// Calculate the position at the entrance of D1
ye = yc - TMath::Abs(theta*r/cs2)*al[1];
xe = xc - r*(Float_t)(TMath::Cos(th1) - TMath::Cos(th2));
// Update track information on positions
track12_p->SetD1bx(xc);
track12_p->SetD1by(yc);
track12_p->SetD1fx(xe);
track12_p->SetD1fy(ye);
track12_p->SetStatus(0);
if(TMath::Abs(xe) > (d1x-t1t2_par[4]) ||
TMath::Abs(ye) > (d1y-t1t2_par[4])) {
track12_p->SetStatus(1);
}
if(TMath::Abs(xc) > (d1x-t1t2_par[4]) ||
TMath::Abs(yc) > (d1y-t1t2_par[4])) {
track12_p->SetStatus(2);
}
// Convert to global coordinates and project to nominal vertex
xl[0] = xe;
xl[1] = ye;
xl[2] = -d1l;
d1->LocalToGlobal(xl,xg,0);
d1->LocalToGlobal(aln,ag,1);
// Project to vertex (z=0)
t = -xg[2] / ag[2];
xt = xg[0] + t * ag[0];
yt = xg[1] + t * ag[1];
track12_p->SetTx(xt);
track12_p->SetTy(yt);
// Calculate scattering angles
xl[0] = xg[0] - xt;
xl[1] = xg[1] - yt;
xl[2] = xg[2];
phi = (Float_t)TMath::ATan(xl[1] / xl[0]) * (Float_t)raddeg;
xt = (Float_t)sqrt(xl[0]*xl[0] + xl[1]*xl[1]);
theta = (Float_t)TMath::ATan(xt / xl[2] ) * (Float_t)raddeg;
track12_p->SetTheta(theta);
track12_p->SetPhi(phi);
}
}
void BrCombineTrack::CombineT3T4T5(BrEventNode *InputTable,BrEventNode *OutputTable)
{
//This routine will independently lookup local tracks in T4 and T3
//and look for combinations which will match in the center of the
//magnet D3. It is required that a track be within the real volume
//of the magnet at this point. In addition a status flag is set if
//in addition is outside a fidicual cut in x and y. The matching has
//cut parameters dx, dy and daly where the position of two tracks must
//be better than in x y and the y direction cosine in order for the
//matching to be accepted.
//
//The same procedure is repeated for T4 and T5. Two sets of new
//objects are created, namely track34 and track45.
//Finally the set of tracks are combined searching for tracks which
//go all the way through the magnet. This is achived by the
//requirement that they have the same T4TR and that they have the
//same momentum to 5% in D3 and D4. These tracks are not complete
//but will be filled with the tracking back to T2 and target
//information.
BrDetectorTrack* t3tr_p;
BrDetectorTrack* t4tr_p;
BrDetectorTrack* t5tr_p;
BrTrack34* track34_p;
BrTrack45* track45_p;
BrTrack345* track345_p;
BrTrack* track_p;
Int_t NumTrack34T4TR;
BrTrack34* track34t4tr[max_t34t4tr];
Float_t d3x,d3y,d3l;
Float_t d4x,d4y,d4l;
Float_t xgt3[3],agt3[3];
Float_t xgt4[3],agt4[3];
Float_t xgt5[3],agt5[3];
static const Int_t max_t3 = 25;
static const Int_t max_t4 = 20;
static const Int_t max_t5 = 15;
Float_t xlt3[max_t3][3],alt3[max_t3][3];
Float_t xlt5[max_t5][3],alt5[max_t5][3];
Float_t xpt3[max_t3],ypt3[max_t3];
Float_t xpt5[max_t5],ypt5[max_t5];
Float_t xpt4d3[max_t4],ypt4d3[max_t4];
Float_t xpt4d4[max_t4],ypt4d4[max_t4];
Float_t xlt4d3[max_t4][3],alt4d3[max_t4][3];
Float_t xlt4d4[max_t4][3],alt4d4[max_t4][3];
Float_t t;
Int_t i,j,it3,it4,it5;
Float_t dx,dy,aly,cs1,cs2,sin1,sin2,p;
BrDetectorVolume *t3;
BrDetectorVolume *t4;
BrDetectorVolume *t5;
BrDetectorVolume *d3;
BrDetectorVolume *d4;
BrDataTable *T3DetectorTracks;
BrDataTable *T4DetectorTracks;
BrDataTable *T5DetectorTracks;
Int_t NumT3Tr,NumT4Tr,NumT5Tr;
Int_t NumTrack45;
T3DetectorTracks = InputTable->GetDataTable("T3 Detector Tracks");
if(T3DetectorTracks == 0) {
cout<<"Can't process CombineT3T4T5 because T3 Detector Tracks table not supplied"<<endl;
}
T4DetectorTracks = InputTable->GetDataTable("T4 Detector Tracks");
if(T4DetectorTracks == 0) {
cout<<"Can't process CombineT3T4T5 because T4 Detector Tracks table not supplied"<<endl;
}
T5DetectorTracks = InputTable->GetDataTable("T5 Detector Tracks");
if(T5DetectorTracks == 0) {
cout<<"Can't process CombineT3T4T5 because T5 Detector Tracks table not supplied"<<endl;
}
if(!T3DetectorTracks || !T4DetectorTracks || !T5DetectorTracks) return;
//Get pointers to volumes we need.
t3 = GetDetectorVolume("T3");
if(!t3) cout<<"Can't process CombineT3T4T5 because T3 volume not supplied"<<endl;
t4 = GetDetectorVolume("T4");
if(!t4) cout<<"Can't process CombineT3T4T5 because T4 volume not supplied"<<endl;
t5 = GetDetectorVolume("T5");
if(!t5) cout<<"Can't process CombineT3T4T5 because T5 volume not supplied"<<endl;
d3 = GetDetectorVolume("D3");
if(!d3) cout<<"Can't process CombineT3T4T5 because D3 volume not supplied"<<endl;
d4 = GetDetectorVolume("D4");
if(!d4) cout<<"Can't process CombineT3T4T5 because D4 volume not supplied"<<endl;
if(!t3 || !t4 || !t5 || !d3 || !d4) return;
NumT3Tr = T3DetectorTracks->GetEntries();
NumT4Tr = T4DetectorTracks->GetEntries();
NumT5Tr = T5DetectorTracks->GetEntries();
//Prepare the output table
Char_t TableName[80];
sprintf(TableName,"BRAHMS Tracksx0");
BrDataTable *Tracks = new BrDataTable(TableName);
OutputTable->AddObject(Tracks);
d3x = d3->GetSize()[0] / (Float_t)2.0;
d3y = d3->GetSize()[1] / (Float_t)2.0;
d3l = d3->GetSize()[2] / (Float_t)2.0;
d4x = d4->GetSize()[0] / (Float_t)2.0;
d4y = d4->GetSize()[1] / (Float_t)2.0;
d4l = d4->GetSize()[2] / (Float_t)2.0;
if(NumT3Tr > max_t3) cout<<"Exceeding max_t3 limit"<<endl;
if(NumT4Tr > max_t4) cout<<"Exceeding max_t4 limit"<<endl;
if(NumT5Tr > max_t5) cout<<"Exceeding max_t5 limit"<<endl;
//Loop over tracks from T4 and get their coordinates
//in D3 and D4 frame
for(i=0;i<NumT4Tr;i++) {
t4tr_p = (BrDetectorTrack*)T4DetectorTracks->At(i);
t4->LocalToGlobal(t4tr_p->GetPos(),xgt4,0);
t4->LocalToGlobal(t4tr_p->GetAlph(),agt4,1);
// Convert into system for magnet D3
d3->GlobalToLocal(xgt4,xlt4d3[i],0);
d3->GlobalToLocal(agt4,alt4d3[i],1);
// Track to center of magnet D3
t = - xlt4d3[i][2] / alt4d3[i][2];
xpt4d3[i] = xlt4d3[i][0] + alt4d3[i][0] * t;
ypt4d3[i] = xlt4d3[i][1] + alt4d3[i][1] * t;
// Convert into system for magnet D4
d4->GlobalToLocal(xgt4,xlt4d4[i],0);
d4->GlobalToLocal(agt4,alt4d4[i],1);
// Track to center of magnet D4
t = - xlt4d4[i][2] / alt4d4[i][2];
xpt4d4[i] = xlt4d4[i][0] + alt4d4[i][0] * t;
ypt4d4[i] = xlt4d4[i][1] + alt4d4[i][1] * t;
}
// Loop over tracks from T3 and track to D3 center
for(i=0;i<NumT3Tr;i++) {
t3tr_p = (BrDetectorTrack*)T3DetectorTracks->At(i);
t3->LocalToGlobal(t3tr_p->GetPos(),xgt3,0);
t3->LocalToGlobal(t3tr_p->GetAlph(),agt3,1);
// Convert into system for magnet D3
d3->GlobalToLocal(xgt3,xlt3[i],0);
d3->GlobalToLocal(agt3,alt3[i],1);
// Track to center of magnet D3
t = - xlt3[i][2] / alt3[i][2];
xpt3[i] = xlt3[i][0] + alt3[i][0] * t;
ypt3[i] = xlt3[i][1] + alt3[i][1] * t;
}
// Loop over tracks from T5 and track to D4 center
for(i=0;i<NumT5Tr;i++) {
t5tr_p = (BrDetectorTrack*)T5DetectorTracks->At(i);;
t5->LocalToGlobal(t5tr_p->GetPos(),xgt5,0);
t5->LocalToGlobal(t5tr_p->GetAlph(),agt5,1);
// Convert into system for magnet D4
d4->GlobalToLocal(xgt5,xlt5[i],0);
d4->GlobalToLocal(agt5,alt5[i],1);
// Track to center of magnet D4
t = - xlt5[i][2] / alt5[i][2];
xpt5[i] = xlt5[i][0] + alt5[i][0] * t;
ypt5[i] = xlt5[i][1] + alt5[i][1] * t;
}
//Track matching part
//Generate T3-T4 matching of tracks
for(it4=0;it4<NumT4Tr;it4++) {
t4tr_p = (BrDetectorTrack*)T4DetectorTracks->At(it4);
if(TMath::Abs(xpt4d3[it4]) < d3x && TMath::Abs(ypt4d3[it4]) < d3y) {
// Accepted in D3, continue looping
for(it3=0;it3<NumT3Tr;it3++) {
t3tr_p = (BrDetectorTrack*)T3DetectorTracks->At(it3);
if(TMath::Abs(xpt3[it3])<d3x && TMath::Abs(ypt3[it3])<d3y) {
dx = xpt4d3[it4] - xpt3[it3];
dy = ypt4d3[it4] - ypt3[it3];
aly = alt4d3[it4][1] - alt3[it3][1];
// If matching condition is fullfilled calculate momentum
if(TMath::Abs(dx) < t345_par[0] &&
TMath::Abs(dy) < t345_par[1] &&
TMath::Abs(aly)< t345_par[2] ) {
// Calculate cos alpha values
cs1 = (Float_t)0.0;
cs2 = (Float_t)0.0;
for(j=0;j<3;j+=2) {
cs1 += alt3[it3][j] * alt3[it3][j];
cs2 += alt4d3[it4][j] * alt4d3[it4][j];
}
cs1 = (Float_t)sqrt(cs1);
cs2 = (Float_t)sqrt(cs2);
sin1 = - alt3[it3][0] / cs1;
sin2 = - alt4d3[it4][0] / cs2;
p = b3 * (d3l + d3l) / ((Float_t)3335.6 * (sin1-sin2));
p = p / (Float_t)sqrt((Float_t)1.0 - alt3[it3][1]*alt3[it3][1]);
BrClonesArray &track34 = *Tracks34;
Int_t ntrack34 = Tracks34->GetEntries();
track34_p = new(track34[ntrack34]) BrTrack34();
track34_p->SetP(p);
track34_p->SetDx(dx);
track34_p->SetDy(dy);
track34_p->SetAly(aly);
track34_p->SetT3tr(t3tr_p);
track34_p->SetT4tr(t4tr_p);
}
}
else {
// Rejected in T3
}
}
}
else {
// Rejected in D3, give up
}
}
//Generate T4-T5 matching of tracks
for(it4=0;it4<NumT4Tr;it4++) {
t4tr_p = (BrDetectorTrack*)T4DetectorTracks->At(it4);
if(TMath::Abs(xpt4d4[it4]) < d4x && TMath::Abs(ypt4d4[it4]) < d4y) {
for(it5=0;it5<NumT5Tr;it5++) {
t5tr_p = (BrDetectorTrack*)T5DetectorTracks->At(it5);
if(TMath::Abs(xpt5[it5])<d4x && TMath::Abs(ypt5[it5])<d4y) {
// This is a possible candidate
dx = xpt4d4[it4] - xpt5[it5];
dy = ypt4d4[it4] - ypt5[it5];
aly = alt4d4[it4][1] - alt5[it5][1];
// If matching condition is fullfilled calculate momentum
if(TMath::Abs(dx) < t345_par[0] &&
TMath::Abs(dy) < t345_par[1] &&
TMath::Abs(aly)< t345_par[2] ) {
// Calculate cos alpha values
cs1 = (Float_t)0.0;
cs2 = (Float_t)0.0;
for(j=0;j<3;j+=2) {
cs1 += alt4d4[it4][j] * alt4d4[it4][j];
cs2 += alt5[it5][j] * alt5[it5][j];
}
cs1 = (Float_t)sqrt(cs1);
cs2 = (Float_t)sqrt(cs2);
sin1 = - alt4d4[it4][0] / cs1;
sin2 = - alt5[it5][0] / cs2;
p = b4 * (d4l + d4l) / ((Float_t)3335.6 * (sin1-sin2));
p = p / (Float_t)sqrt((Float_t)1.0 - alt5[it5][1]*alt5[it5][1]);
BrClonesArray &track45 = *Tracks45;
Int_t ntrack45 = Tracks45->GetEntries();
track45_p = new(track45[ntrack45]) BrTrack45();
track45_p->SetP(p);
track45_p->SetDx(dx);
track45_p->SetDy(dy);
track45_p->SetAly(aly);
track45_p->SetT4tr(t4tr_p);
track45_p->SetT5tr(t5tr_p);
}
}
else {
// Rejected in D4
}
}
}
else {
// Rejected in D4, give up
}
}
//Combine track34 and track45
NumTrack45 = Tracks45->GetEntries();
for(it4=0;it4<NumTrack45;it4++) {
track45_p = (BrTrack45*)Tracks45->At(it4);
t4tr_p = track45_p->GetT4tr();
// Need list of track34 associated with track45
// ie look for track34 that have same t4tr as track45
GetTrack34WithSameT4TR(t4tr_p,track34t4tr,&NumTrack34T4TR);
for(it3=0;it3<NumTrack34T4TR;it3++) {
track34_p = track34t4tr[it3];
if(TMath::Abs((track34_p->GetP()-track45_p->GetP())/track34_p->GetP()) < .10) {
BrClonesArray &track345 = *Tracks345;
Int_t ntrack345 = Tracks345->GetEntries();
track345_p = new(track345[ntrack345]) BrTrack345();
ntrack345 = Tracks345->GetEntries();
track345_p->SetP((track34_p->GetP() + track45_p->GetP()) / (Float_t)2.0);
track34_p->SetTrack345(track345_p);
track45_p->SetTrack345(track345_p);
track345_p->SetTrack34(track34_p);
track345_p->SetTrack45(track45_p);
// Create the first information for the BRAHMS track
track_p = new BrTrack();
Tracks->Add(track_p);
track_p->SetP34(track34_p->GetP());
track_p->SetP45(track45_p->GetP());
track_p->SetP345(track345_p->GetP());
track34_p->SetTrack(track_p);
track45_p->SetTrack(track_p);
track345_p->SetTrack(track_p);
}
}
}
}
void BrCombineTrack::GetTrack34WithSameT4TR(BrDetectorTrack* t4tr_p,BrTrack34 *strack[max_t34t4tr],Int_t* NumFound)
{
//Loop over tracks that were combined in T3 and T4 and pick the ones that have the
//same Detector Track in T4 as in the argument. Routine meant to be used to look
//for track34's that have same t4tr as track45.
Int_t i,numfound;
Int_t NumTrack34;
BrTrack34 *track34_p;
numfound = 0;
NumTrack34 = Tracks34->GetEntries();
for(i=0;i<NumTrack34;i++) {
track34_p = (BrTrack34*)Tracks34->At(i);
if( track34_p->GetT4tr() == t4tr_p ) {
strack[numfound] = track34_p;
numfound++;
}
}
*NumFound = numfound;
}
void BrCombineTrack::CombineD12D34(BrEventNode *InputTable,BrEventNode *OutputTable)
{
//This routine will take tracks from T3/T4/T5 and project them
//to T2 and find candidates for the track.
//If found it will update parts of the real track information.
BrTrack12* track12_p;
BrTrack34* track34_p;
BrTrack345* track345_p;
BrTrack* track_p;
BrDetectorTrack* t2tr_p;
BrDetectorTrack* t3tr_p;
BrT2Match* t2match_p;
Float_t xg[3],ag[3];
Float_t xl[3],al[3];
Float_t t,xt2,yt2,dx,dy,daly;
Int_t it,it2;
BrDetectorVolume *t2;
BrDetectorVolume *t3;
BrDataTable *T1DetectorTracks;
BrDataTable *T2DetectorTracks;
BrDataTable *Tracks;
Int_t NumT2Tr,NumTracks345;
t2 = GetDetectorVolume("T2");
if(!t2) cout<<"Can't process CombineD12D34 because T2 volume not supplied"<<endl;
t3 = GetDetectorVolume("T3");
if(!t3) cout<<"Can't process CombineD12D34 because T3 volume not supplied"<<endl;
T1DetectorTracks = InputTable->GetDataTable("T1 Detector Tracks");
if(!T1DetectorTracks) cout<<"Can't process CombineD12D34 because T1 Detector Tracks not supplied"<<endl;
T2DetectorTracks = InputTable->GetDataTable("T2 Detector Tracks");
if(!T2DetectorTracks) cout<<"Can't process CombineD12D34 because T2 Detector Tracks not supplied"<<endl;
Tracks = OutputTable->GetDataTable("BRAHMS Tracks");
if(!Tracks) cout<<"Can't process CombineD12D34 because no BRAHMS Tracks produced"<<endl;
if(!T1DetectorTracks||!T2DetectorTracks||!Tracks) return;
NumT2Tr = T2DetectorTracks->GetEntries();
NumTracks345 = Tracks345->GetEntries();
//Loop over T345 track entries
for(it=0;it<NumTracks345;it++) {
track345_p = (BrTrack345*)Tracks345->At(it);;
track34_p = track345_p->GetTrack34();
track_p = track345_p->GetTrack(); //Get the track that we started to construct
t3tr_p = track34_p->GetT3tr();
// Get track coordinates in T3 and transform to T2 frame
t3->LocalToGlobal(t3tr_p->GetPos(),xg,0);
t3->LocalToGlobal(t3tr_p->GetAlph(),ag,1);
t2->GlobalToLocal(xg,xl,0);
t2->GlobalToLocal(ag,al,1);
// Track to center of T2
t = - xl[2] / al[2];
xt2 = xl[0] + t * al[0];
yt2 = xl[1] + t * al[1];
track345_p->SetT2x(xt2);
track345_p->SetT2y(yt2);
// Loop over possible T2TR;
for(it2=0;it2<NumT2Tr;it2++) {
t2tr_p = (BrDetectorTrack*)T2DetectorTracks->At(it2);
track12_p = GetTrack12WithSameT2Tr(t2tr_p);
dx = t2tr_p->GetPos()[0] - xt2;
dy = t2tr_p->GetPos()[1] - yt2;
daly = t2tr_p->GetAlph()[1] - al[1];
// Set criteria for matching to tracks
// one could select to let the T2TR tracks point to the
// real track and then do the selection
if(TMath::Abs(dx) < t2m_par[0] &&
TMath::Abs(dy) < t2m_par[1] &&
TMath::Abs(daly) < t2m_par[2] ) {
if(!track12_p) {
cout<<"track12_p == 0 and it should not!!!!!!!!!!!!!"<<endl;
cout<<"Expect this to crash very soon!!!!!!"<<endl;
}
BrClonesArray &t2match = *T2Match;
Int_t nt2match = T2Match->GetEntries();
t2match_p = new(t2match[nt2match]) BrT2Match();
t2match_p->SetDelx(dx);
t2match_p->SetDely(dy);
t2match_p->SetDelay(daly);
t2match_p->SetTx((Float_t)0.0);
t2match_p->SetTy((Float_t)0.0);
t2match_p->SetD1fx((Float_t)0.0);
t2match_p->SetD1fy((Float_t)0.0);
t2match_p->SetD2fx((Float_t)0.0);
t2match_p->SetD2fy((Float_t)0.0);
t2match_p->SetQval((Float_t)sqrt(dx*dx + dy*dy+(Float_t)10000.*daly*daly));
t2match_p->SetT2tr(t2tr_p);
t2match_p->SetTrack12(track12_p);
t2match_p->SetTrack345(track345_p);
t2match_p->SetTrack(track_p);
}
}
}
}
BrTrack12 *BrCombineTrack::GetTrack12WithSameT2Tr(BrDetectorTrack *t2tr_p){
for(Int_t i=0;i<Tracks12->GetEntries();i++) {
BrTrack12 *track12_p = (BrTrack12*)Tracks12->At(i);
if(track12_p->GetT2tr() == t2tr_p) return track12_p;
}
return 0;
}
void BrCombineTrack::Track345ToTarget()
{
//This routine will loop over all tracks found using their position
//in T2 as a starting point for backward tracking. It will
//the Tx and Ty at the target position, calculate the scattering
//angle using the information for the track positions at the
//entrance to D1, ie on a virtual plane. The momentum is taken
//from the values obtained from the T3/T4/T5 tracking.
BrDetectorTrack* t2tr_p;
BrTrack12 *track12_p;
BrTrack345 *track345_p;
BrTrack *track_p;
BrT2Match *t2match_p;
Float_t d1x,d1y,d1l;
Float_t d2x,d2y,d2l;
Float_t xg[3],ag[3];
Float_t xl[3],al[3];
Float_t aln[3];
Float_t t,xc,yc,cs2,r;
Float_t sin1,sin2,th1,th2;
Float_t theta,phi,cs,si,xe,ye,xd1,yd1,xt,yt;
Int_t i,it;
Int_t NumT2Match;
BrDetectorVolume *t1;
BrDetectorVolume *t2;
BrDetectorVolume *d1;
BrDetectorVolume *d2;
t1 = GetDetectorVolume("T1");
if(!t1) cout<<"Can't process Track345ToTarget because T1 volume not supplied"<<endl;
t2 = GetDetectorVolume("T2");
if(!t2) cout<<"Can't process Track345ToTarget because T2 volume not supplied"<<endl;
d1 = GetDetectorVolume("D1");
if(!d1) cout<<"Can't process Track345ToTarget because D1 volume not supplied"<<endl;
d2 = GetDetectorVolume("D2");
if(!d2) cout<<"Can't process Track345ToTarget because D2 volume not supplied"<<endl;
if(!t1 || !t2 || !d1 || !d2) return;
d1x = t1->GetSize()[0] / (Float_t)2.0;
d1y = t1->GetSize()[1] / (Float_t)2.0;
d1l = t1->GetSize()[2] / (Float_t)2.0;
d2x = t2->GetSize()[0] / (Float_t)2.0;
d2y = t2->GetSize()[1] / (Float_t)2.0;
d2l = t2->GetSize()[2] / (Float_t)2.0;
//Loop over tracks
NumT2Match = T2Match->GetEntries();
for(it=0;it<NumT2Match;it++) {
t2match_p = (BrT2Match*)T2Match->At(it);
t2tr_p = t2match_p->GetT2tr();
track12_p = t2match_p->GetTrack12();
track345_p = t2match_p->GetTrack345();
track_p = t2match_p->GetTrack();
t2->LocalToGlobal(t2tr_p->GetPos(),xg,0);
t2->LocalToGlobal(t2tr_p->GetAlph(),ag,1);
d2->GlobalToLocal(xg,xl,0);
d2->GlobalToLocal(ag,al,1);
// Project to front of D2 and bend according to value of p
// and magnetic fields.
t = -(xl[2] - d2l) / al[2];
xc = xl[0] + t * al[0];
yc = xl[1] + t * al[1];
cs2 = (Float_t)sqrt(al[0]*al[0] + al[2]*al[2]);
sin2 = al[0] / cs2;
r = (Float_t) 3335.6 * track345_p->GetP() * cs2 / b2;
sin1 = -(Float_t)2.0 * d2l / r + sin2;
th1 = (Float_t)TMath::ASin(sin1);
th2 = (Float_t)TMath::ASin(sin2);
theta = th2 - th1;
cs = (Float_t)TMath::Cos(theta);
si = (Float_t)TMath::Sin(theta);
// Rotate the coordinates for (al[1],al[2]) by the scattering angle
aln[0] = al[0] * cs - al[2] * si;
aln[2] = al[0] * si + al[2] * cs;
aln[1] = al[1];
// These are now direction cosines.
// Calculate the position at the entrance of D2
ye = yc - TMath::Abs(theta * r / cs2) * al[1];
xe = xc + r*(Float_t)(TMath::Cos(th2) - TMath::Cos(th1));
xl[0] = xe;
xl[1] = ye;
xl[2] = -d2l;
t2match_p->SetD2bx(xc);
t2match_p->SetD2by(yc);
t2match_p->SetD2fx(xe);
t2match_p->SetD2fy(ye);
// We now have coordinates in D2 system, go to D1
d2->LocalToGlobal(xl,xg,0);
d2->LocalToGlobal(aln,ag,1);
// Project to center of D1 or front depending on algorithm
t = - (xl[2] - d1l) / al[2];
xc = xl[0] + t * al[0];
yc = xl[1] + t * al[1];
cs2 = (Float_t)sqrt(al[0]*al[0] + al[2]*al[2]);
sin2 = al[0] / cs2;
r = (Float_t)3335.6 * track345_p->GetP() * cs2 / b1;
sin1 = -(Float_t)2.0 * d1l/r + sin2;
th1 = (Float_t)TMath::ASin(sin1);
th2 = (Float_t)TMath::ASin(sin2);
theta = th2 - th1;
// Rotate the coordinates for al[0],al[2] by the scattering angle
cs = (Float_t)TMath::Cos(theta);
si = (Float_t)TMath::Sin(theta);
aln[0] = al[0] * cs - al[2] * si;
aln[2] = al[0] * si + al[2] * cs;
aln[1] = al[1];
// These are now the direction cosines. Calculate the position
// at the entrance of D1
yd1 = yc - TMath::Abs(theta*r/cs2) * al[1];
xd1 = xc + r * (Float_t)(TMath::Cos(th2) - TMath::Cos(th1));
t2match_p->SetD1bx(xc);
t2match_p->SetD1by(yc);
t2match_p->SetD1fx(xd1);
t2match_p->SetD1fy(yd1);
// Convert center point to global coordinates
xl[0] = xd1;
xl[1] = yd1;
xl[2] = -d1l;
d1->LocalToGlobal(xl,xg,0);
d1->LocalToGlobal(aln,ag,1);
// Project to target (z=0)
t = -xg[2] / ag[2];
xt = xg[0] + t * ag[0];
yt = xg[1] + t * ag[1];
t2match_p->SetTx(xt);
t2match_p->SetTy(yt);
// Calculate scattering angles
for(i=0;i<3;i++) xl[i] = xg[i];
phi = (Float_t)TMath::ATan(xl[1] / xl[0]) * (Float_t)raddeg;
xt = (Float_t)sqrt(xl[0]*xl[0] + xl[1]*xl[1]);
theta = (Float_t)TMath::ATan(xt / xl[2] ) * (Float_t)raddeg;
t2match_p->SetTheta(theta);
t2match_p->SetPhi(phi);
// Now put all relevant information into track that will be returned.
track_p->SetP12(track12_p->GetP());
track_p->SetP((track12_p->GetP() + (Float_t)2.0 * track_p->GetP345())/(Float_t)3.0);
track_p->SetTheta(theta);
track_p->SetPhi(phi);
}
}
void BrCombineTrack::CombineMIDST1T2(BrEventNode *InputTable)
{
//Loop over tracks in T1/T2 and see how well they match and
//create candidates for tracks
//This is a first primitive attempt on match tracks through
//the magnet M0.
BrDetectorTrack* t1tr_p;
BrDetectorTrack* t2tr_p;
BrTrack12* track12_p;
Int_t i,j;
Float_t m0x,m0y,m0l,t;
Float_t xgt1[3],agt1[3];
Float_t xgt2[3],agt2[3];
#define max_t1 100
#define max_t2 50
Float_t xlt1[max_t1][3],alt1[max_t1][3];
Float_t xlt2[max_t2][3],alt2[max_t2][3];
Float_t xpt1[max_t1],ypt1[max_t1],zpt1[max_t1];
Float_t xpt2[max_t2],ypt2[max_t2],zpt2[max_t2];
Float_t xct1[max_t1],yct1[max_t1];
Float_t xct2[max_t1],yct2[max_t1];
Float_t cs1,cs2,sin1,sin2,dx,dy,aly,p,p_xz;
Float_t theta_enter,theta_exit,theta_link,ang_enter,ang_exit;
Float_t dang;
Int_t status_t1t2;
Int_t it1,it2;
BrDetectorVolume *m0;
BrDetectorVolume *t1;
BrDetectorVolume *t2;
BrDataTable *T1DetectorTracks;
BrDataTable *T2DetectorTracks;
Int_t NumT1Tr,NumT2Tr;
T1DetectorTracks = InputTable->GetDataTable("T1 Detector Tracks");
if(!T1DetectorTracks) cout<<"Can't process CombineMIDST1T2 because T1 Detector Tracks volume not supplied"<<endl;
T2DetectorTracks = InputTable->GetDataTable("T2 Detector Tracks");
if(!T2DetectorTracks) cout<<"Can't process CombineMIDST1T2 because T2 Detector Tracks volume not supplied"<<endl;
if(!T1DetectorTracks || !T2DetectorTracks) return;
m0 = GetDetectorVolume("M0");
if(!m0) cout<<"Can't process CombineMIDST1T2 because M0 volume not supplied"<<endl;
t1 = GetDetectorVolume("T1");
if(!t1) cout<<"Can't process CombineMIDST1T2 because T1 volume not supplied"<<endl;
t2 = GetDetectorVolume("T2");
if(!t2) cout<<"Can't process CombineMIDST1T2 because T2 volume not supplied"<<endl;
if(!m0 || !t1 || !t2) return;
//This simple version will loop over the tracks found in T2 and
//convert the local to global coordinates, project the track to
//the center of the magnet.
//It will project the tracks in T1 forward to the center plane of
//the magnet. It will do the calculation in the local coordinate
//system of D2 which is most useful for this purpose.
//It will make simple corrections due to path length differences.
//Find the M0 sizes in the first time
m0x = m0->GetSize()[0] / (Float_t)2.0;
m0y = m0->GetSize()[1] / (Float_t)2.0;
m0l = m0->GetSize()[2] / (Float_t)2.0;
//Loop over tracks from t2
NumT2Tr = T2DetectorTracks->GetEntries();
for(i=0;i<NumT2Tr;i++) {
t2tr_p = (BrDetectorTrack*)T2DetectorTracks->At(i);;
t2->LocalToGlobal(t2tr_p->GetPos(),xgt2,0);
t2->LocalToGlobal(t2tr_p->GetAlph(),agt2,1);
// Convert into system for magnet M0
m0->GlobalToLocal(xgt2,xlt2[i],0);
m0->GlobalToLocal(agt2,alt2[i],1);
// Track to center of effective edge of magnet M0
t = (xlt2[i][2] - m0l )/ alt2[i][2];
xpt2[i] = xlt2[i][0] - alt2[i][0] * t;
ypt2[i] = xlt2[i][1] - alt2[i][1] * t;
zpt2[i] = m0l;
t = -xlt2[i][2] / alt2[i][2];
xct2[i] = xlt2[i][0] + alt2[i][0] * t;
yct2[i] = xlt2[i][1] + alt2[i][1] * t;
}
//Loop over tracks from t1
NumT1Tr = T1DetectorTracks->GetEntries();
for(i=0;i<NumT1Tr;i++) {
t1tr_p = (BrDetectorTrack*)T1DetectorTracks->At(i);
t1->LocalToGlobal(t1tr_p->GetPos(),xgt1,0);
t1->LocalToGlobal(t1tr_p->GetAlph(),agt1,1);
// Convert into system for magnet M0
m0->GlobalToLocal(xgt1,xlt1[i],0);
m0->GlobalToLocal(agt1,alt1[i],1);
// Track to center of effective edge of M0
t = -( m0l + xlt1[i][2] )/ alt1[i][2];
xpt1[i] = xlt1[i][0] + alt1[i][0] * t;
ypt1[i] = xlt1[i][1] + alt1[i][1] * t;
zpt1[i] = -m0l;
t = -xlt1[i][2] / alt1[i][2];
xct1[i] = xlt1[i][0] + alt1[i][0] * t;
yct1[i] = xlt1[i][1] + alt1[i][1] * t;
}
//Loop over combinations of T1 and T2 tracks and look at their
//matching. This is done using the effective edge approximation
//matching on the plane perpendicular to the line connecting the
//two entrance points. The information known is that the two points
//xpt1,xpt2 and the direction cosines. The actual matching is done
//in two subroutines ala E866.
for(it1=0;it1<NumT1Tr;it1++) {
t1tr_p = (BrDetectorTrack*)T1DetectorTracks->At(it1);
// Does the track project into the GAP?
if(TMath::Abs(xpt1[it1])<m0x&&TMath::Abs(ypt1[it1])<m0y) {
// Loop over T2 tracks
for(it2=0;it2<NumT2Tr;it2++) {
t2tr_p = (BrDetectorTrack*)T2DetectorTracks->At(it2);
if(TMath::Abs(xpt2[it2])<m0x&&TMath::Abs(ypt2[it2])<m0y) {
// Calculate the approximate momentum
// (to be used for dx,dy cuts)
cs1 = (Float_t)0.0;
cs2 = (Float_t)0.0;
for(j=0;j<3;j+=2) {
cs1 += alt1[it1][j] * alt1[it1][j];
cs2 += alt2[it2][j] * alt2[it2][j];
}
cs1 = (Float_t)sqrt(cs1);
cs2 = (Float_t)sqrt(cs2);
sin1 = alt1[it1][0] / cs1;
sin2 = alt2[it2][0] / cs2;
p_xz = b0 * (m0l+m0l) / ((Float_t)3335.6*(sin1-sin2));
theta_enter = (Float_t)TMath::ASin(sin1);
theta_exit = (Float_t)TMath::ASin(sin2);
// Obtain angles used for cuts
theta_link = (Float_t)TMath::ATan((xpt2[it2]-xpt1[it1])/(2.0*m0l));
ang_enter = theta_enter - theta_link;
ang_exit = theta_exit - theta_link;
dang = ang_enter + ang_exit;
dy = yct2[it2] - yct1[it1];
aly = alt2[it2][1] - alt1[it1][1];
// Match tracks in theta1-theta2 and then perform the
// match in the other plane
// Match tracks and calculate momentum
if(TMath::Abs(dang) < t1t2_par[3] &&
TMath::Abs(dy) < t1t2_par[1] &&
TMath::Abs(aly) < t1t2_par[3] ) {
// Set some status flag
status_t1t2 = 0;
if((TMath::Abs(xpt2[it2]) > (m0x-t1t2_par[0]))||
(TMath::Abs(ypt2[it2]) > (m0y-t1t2_par[1]))) {
status_t1t2 |= 1<<sf_t1_d2_midcut;
}
if((TMath::Abs(xpt1[it1]) > (m0x-t1t2_par[0]))||
(TMath::Abs(ypt1[it1]) > (m0y-t1t2_par[1]))) {
status_t1t2 |= 1<<sf_t2_d2_midcut;
}
p = p_xz / (Float_t)TMath::Sqrt(1.0-alt1[it1][1]*alt1[it1][1]);
BrClonesArray &track12 = *Tracks12;
Int_t ntrack12 = Tracks12->GetEntries();
track12_p = new(track12[ntrack12]) BrTrack12();
track12_p->SetP(p);
track12_p->SetDx(dx);
track12_p->SetDy(dy);
track12_p->SetT1tr(t1tr_p);
track12_p->SetT2tr(t2tr_p);
}
}
}
}
else {
// t1 track rejected, put appropriate stuff here
}
}
}
void BrCombineTrack::Track12MIDSToTarget()
{
//This routine will loop over all tracks found using their position
//in T1 as a starting point for backward tracking. It will
//calculate Tx and Ty at the target position, calculate the
//scattering angle using the information for the track positions
//TPC1.
//The target Tx and Ty are in the SPEC coordinate system ie at
//z' = 0 (ie. on the beam axis parallel to the M0 front face)
BrDetectorTrack* t1tr_p;
BrTrack12* track12_p;
Int_t i,it;
Float_t x0[3],xl[3],al[3],xg[3],ag[3];
Float_t t,xt,yt,theta,phi;
Int_t NumTrack12;
BrDetectorVolume *m0;
BrDetectorVolume *t1;
m0 = GetDetectorVolume("M0");
if(!m0) cout<<"Can't process Track12MIDSToTarget because M0 volume not supplied"<<endl;
t1 = GetDetectorVolume("T1");
if(!t1) cout<<"Can't process Track12MIDSToTarget because T1 volume not supplied"<<endl;
//Find the origin of spectrometer in SPEC system
for(i=0;i<3;i++) xg[i] = (Float_t)0.0;
m0->GlobalToLocal(xg,x0,0);
//Loop over tracks
NumTrack12 = Tracks12->GetEntries();
for(it=0;it<NumTrack12;it++) {
track12_p = (BrTrack12*)Tracks12->At(it);
t1tr_p = track12_p->GetT1tr();
// The local track coordinates are given in the local coordinates
// for T1. Convert to global --> local for D1
t1->LocalToGlobal(t1tr_p->GetPos(),xg,0);
t1->LocalToGlobal(t1tr_p->GetAlph(),ag,1);
m0->GlobalToLocal(xg,xl,0);
m0->GlobalToLocal(ag,al,1);
// Update track information on positions
track12_p->SetD1bx((Float_t)0.0);
track12_p->SetD1by((Float_t)0.0);
track12_p->SetD1fx((Float_t)0.0);
track12_p->SetD1fy((Float_t)0.0);
track12_p->SetStatus(0);
// Project to vertex (z=0)
t = (x0[2] - xl[2] ) / al[2];
xt = xl[0] + t * al[0];
yt = xl[1] + t * al[1];
track12_p->SetTx(xt);
track12_p->SetTy(yt);
// Calculate scattering angles
xl[0] = xl[0] - xt;
xl[1] = xl[1] - yt;
xl[2] = xl[2] - x0[2];
m0->LocalToGlobal(xl,xg,1);
phi = (Float_t)TMath::ATan(xg[1] / xg[0]) * (Float_t)raddeg;
xt = (Float_t)sqrt(xg[0]*xg[0] + xg[1]*xg[1]);
theta = (Float_t)TMath::ATan(xt / xg[2] ) * (Float_t)raddeg;
track12_p->SetTheta(theta);
track12_p->SetPhi(phi);
}
}
Int_t BrCombineTrack::GetTrackingLevel(Char_t *detector)
{
return kTRUE; //do that for now, make more sophisticated later
}
void BrCombineTrack::ReadSimulationParameters( Char_t *filename )
{
FILE *ParFile;
#define MaxLine 80
Char_t label[MaxLine];
ParFile = fopen( filename,"r");
if(!ParFile) {
cout<<"Error opening Simulation Parameter File "<<filename<<endl;
cout<<"Please set the filename correctly!!!"<<endl;
cout<<"Returning from this routine without doing anything!!!"<<endl;
fSimulationParametersRead = kFALSE;
return;
}
fgets(label,MaxLine,ParFile);
fgets(label,MaxLine,ParFile);
sscanf(label,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d",
&pr_opt[0],&pr_opt[1],&pr_opt[2],&pr_opt[3],&pr_opt[4],&pr_opt[5],&pr_opt[6],
&pr_opt[7],&pr_opt[8],&pr_opt[9],&pr_opt[10],&pr_opt[11],&pr_opt[12],&pr_opt[13],
&pr_opt[14],&pr_opt[15],&pr_opt[16],&pr_opt[17],&pr_opt[18],&pr_opt[19]);
fgets(label,MaxLine,ParFile);
fgets(label,MaxLine,ParFile);
sscanf(label,"%f %f %f %f %f %f %f %f %f %f",
&t1_par[0],&t1_par[1],&t1_par[2],&t1_par[3],&t1_par[4],
&t1_par[5],&t1_par[6],&t1_par[7],&t1_par[8],&t1_par[9]);
fgets(label,MaxLine,ParFile);
fgets(label,MaxLine,ParFile);
sscanf(label,"%f %f %f %f %f %f %f %f %f %f",
&t2_par[0],&t2_par[1],&t2_par[2],&t2_par[3],&t2_par[4],
&t2_par[5],&t2_par[6],&t2_par[7],&t2_par[8],&t2_par[9]);
fgets(label,MaxLine,ParFile);
fgets(label,MaxLine,ParFile);
sscanf(label,"%f %f %f %f %f %f %f %f %f %f",
&h1_par[0],&h1_par[1],&h1_par[2],&h1_par[3],&h1_par[4],
&h1_par[5],&h1_par[6],&h1_par[7],&h1_par[8],&h1_par[9]);
fgets(label,MaxLine,ParFile);
fgets(label,MaxLine,ParFile);
sscanf(label,"%f %f %f %f %f %f %f %f %f %f",
&t1t2_par[0],&t1t2_par[1],&t1t2_par[2],&t1t2_par[3],&t1t2_par[4],
&t1t2_par[5],&t1t2_par[6],&t1t2_par[7],&t1t2_par[8],&t1t2_par[9]);
fgets(label,MaxLine,ParFile);
fgets(label,MaxLine,ParFile);
sscanf(label,"%f %f %f %f %f %f %f %f %f %f",
&t3_par[0],&t3_par[1],&t3_par[2],&t3_par[3],&t3_par[4],
&t3_par[5],&t3_par[6],&t3_par[7],&t3_par[8],&t3_par[9]);
fgets(label,MaxLine,ParFile);
fgets(label,MaxLine,ParFile);
sscanf(label,"%f %f %f %f %f %f %f %f %f %f",
&t4_par[0],&t4_par[1],&t4_par[2],&t4_par[3],&t4_par[4],
&t4_par[5],&t4_par[6],&t4_par[7],&t4_par[8],&t4_par[9]);
fgets(label,MaxLine,ParFile);
fgets(label,MaxLine,ParFile);
sscanf(label,"%f %f %f %f %f %f %f %f %f %f",
&t5_par[0],&t5_par[1],&t5_par[2],&t5_par[3],&t5_par[4],
&t5_par[5],&t5_par[6],&t5_par[7],&t5_par[8],&t5_par[9]);
fgets(label,MaxLine,ParFile);
fgets(label,MaxLine,ParFile);
sscanf(label,"%f %f %f %f %f %f %f %f %f %f",
&t345_par[0],&t345_par[1],&t345_par[2],&t345_par[3],&t345_par[4],
&t345_par[5],&t345_par[6],&t345_par[7],&t345_par[8],&t345_par[9]);
fgets(label,MaxLine,ParFile);
fgets(label,MaxLine,ParFile);
sscanf(label,"%f %f %f %f %f %f %f %f %f %f",
&t2m_par[0],&t2m_par[1],&t2m_par[2],&t2m_par[3],&t2m_par[4],
&t2m_par[5],&t2m_par[6],&t2m_par[7],&t2m_par[8],&t2m_par[9]);
fSimulationParametersRead = kTRUE;
fclose(ParFile);
}
void BrCombineTrack::ReadFieldParameters( Char_t *filename )
{
FILE* FieldFile = fopen( filename,"r");
Char_t line[80];
fgets(line,80,FieldFile);
sscanf(line,"%f %f %f %f",&b1,&b2,&b3,&b4);
fgets(line,80,FieldFile);
sscanf(line,"%f",&b0);
fclose(FieldFile);
}
void BrCombineTrack::AddDetectorVolume(BrDetectorVolume* vol)
{
if(!vol) return;
if(NumDetectorVolumes > max_num_volumes) {
cout<<"Too many volumes specified, need to increase max_num_volumes and recompile"<<endl;
return;
}
fVolume_p[NumDetectorVolumes++] = vol;
}
BrDetectorVolume* BrCombineTrack::GetDetectorVolume(const Char_t *vol) const
{
Int_t i;
for(i=0;i<NumDetectorVolumes;i++) {
if(!strcasecmp(fVolume_p[i]->GetName(),vol)) return fVolume_p[i];
}
return 0;
}
void BrCombineTrack::ListDetectorParameters()
{
cout<<"This routine needs to be implemented"<<endl;
}
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.