mTecPIDModule.cc
#include <stdio.h>
#include <math.h>
#include <algo.h>
#include <vector.h>
#include <map.h>
#include <stdlib.h>
#include <iostream.h>
#include "gsl/gsl_rng.h"
#include "PHNode.h"
#include "PHIODataNode.h"
#include "PHNodeIterator.h"
#include "mTecPIDModule.h"
#include "TecOutV1.hh"
#include "TecOutV2.hh"
#include "TecTrack.hh"
#include "dPadClusterWrapper.h"
#include "dDchTracksWrapper.h"
#include "dCglTrackWrapper.h"
#include "dCglParticleWrapper.h"
#include "dTofReconstructedWrapper.h"
#include "fkinWrapper.h"
#include "tecghitWrapper.h"
#include "EventHeader.h"
#include "RunHeader.h"
#include "TecCalibrationObject.hh"
#include "DchTrack.h"
#include "PadCluster.h"
#include "CglTrack.h"
#include "PHTrackOut.h"
#include "BbcOut.h"
#include "ZdcOut.h"
#include "VtxOut.h"
#include "TofOut.h"
#include "CrkRing.h"
#include "EmcClusterLocalExt.h"
#include "utiCentrality.h"
#include "mTecUtilities.h"
#include "TecCalibrationObject.hh"
#include "PhHistogramFactory.hh"
#include <TROOT.h>
#include <TFile.h>
#include <TH1.h>
#include <TNtuple.h>
typedef PHIODataNode<dPadClusterWrapper> dPadClusterNode_t;
typedef PHIODataNode<dDchTracksWrapper> dDchTracksNode_t;
typedef PHIODataNode<dCglTrackWrapper> dCglTrackNode_t;
typedef PHIODataNode<dCglParticleWrapper> dCglParticleNode_t;
typedef PHIODataNode<fkinWrapper> fkinNode_t;
typedef PHIODataNode<tecghitWrapper> tecghitNode_t;
typedef PHIODataNode < TObject > ObjectNode_t;
typedef PHIODataNode < TecOut > TecOutNodeShort_t;
typedef PHIODataNode<TecOutV1> TecOutNode_t;
// list of geant particles hitting tec
struct tracklist {
int true_track;
float xinglo[10];
float yinglo[10];
float zinglo[10];
int nghits;
} gtrklist[1000];
static gsl_rng *rng;
static int init_done = 0;
static int evtnum = 0;
//===========================================================================
int mTecPIDModule::GetBeg12(TecGeometryObject* TGO,
float* beg1, float* beg2) {
// Find coordinates of the midddle wire of the first and
// last planes in the south side of sector 1.
// This is temporary. Eventually we will calculate track
// length using intersections of the track with first and last
// planes in corresponding sector.
*beg1 = TGO->getGlobalX(12,TecGeometryObject::get_NumWires(0)/2);
// use plane 3 for "out" coordinate 03/06/02
*beg2 = TGO->getGlobalX(18,TecGeometryObject::get_NumWires(3)/2);
return 0;
}
//===========================================================================
int mTecPIDModule::MakeGTrackList(PHCompositeNode *root) {
int firstplane=0; int lastplane=3;
PHNodeIterator i(root);
//------ Get Geant Tables ---------------------------------------------
fkinNode_t* fkN = static_cast<fkinNode_t*>(i.findFirst("PHIODataNode", "fkin"));
if (!fkN) {
cerr << "mTecPIDModule -> ERROR: fkin table not found." << endl;
return False;
}
tecghitNode_t* tgN = static_cast<tecghitNode_t*>(i.findFirst("PHIODataNode", "tecghit"));
if (!tgN) {
cerr << "mTecPIDModule -> ERROR: tecghit table not found." << endl;
return False;
}
tecghitWrapper* tecghit = tgN->getData();
int gtrknum,trk,was,ighit,j;
gtrknum=0; // number of good geant tracks
for(j=0; j<1000; j++) gtrklist[j].nghits=0;
for (ighit=0; ighit<(int)tecghit->RowCount(); ighit++) {
if(tecghit->get_plane(ighit)>=firstplane &&
tecghit->get_plane(ighit)<=lastplane) {
trk=tecghit->get_mctrack(ighit);
was=-1;
for(j=0; j < gtrknum; j++)
if(trk==gtrklist[j].true_track) {was=j; break;}
if(was==-1) {
gtrklist[gtrknum].true_track = trk;
gtrklist[gtrknum].xinglo[gtrklist[gtrknum].nghits]=tecghit->get_xyzinglo(0,ighit);
gtrklist[gtrknum].yinglo[gtrklist[gtrknum].nghits]=tecghit->get_xyzinglo(1,ighit);
gtrklist[gtrknum].zinglo[gtrklist[gtrknum].nghits]=tecghit->get_xyzinglo(2,ighit);
gtrklist[gtrknum].nghits++;
if(gtrklist[gtrknum].nghits>9) gtrklist[gtrknum].nghits=9;
gtrknum++;
if(gtrknum>999) gtrknum=999;
}
else {
gtrklist[was].xinglo[gtrklist[was].nghits]=tecghit->get_xyzinglo(0,ighit);
gtrklist[was].yinglo[gtrklist[was].nghits]=tecghit->get_xyzinglo(1,ighit);
gtrklist[was].zinglo[gtrklist[was].nghits]=tecghit->get_xyzinglo(2,ighit);
gtrklist[was].nghits++;
if(gtrklist[gtrknum].nghits>9) gtrklist[gtrknum].nghits=9;
}
}
} // ighit loop over tecghit entries
return gtrknum;
}
int mTecPIDModule::GetBetaGamma(float chanperbin, float dedx,
float *betagamma_ptr, int *nbtgm) {
// Calculates beta*gamma from dE/dX assuming that the particle is not
// an electron. Two beta*gamma values are possible.
int i,nbg;
float x,tmp[1000],bg[2];
nbg=0; bg[0]=0.; bg[1]=0.;
for(i=0; i<1000; i++)
{
x=0.1*(i+1);
tmp[i] = chanperbin*0.0705*(1+1/x/x)*(11.45+log(x*x)-(x*x)/(1+x*x));
}
for(i=0; i<999; i++)
{
if(((dedx<tmp[i]) && (dedx>tmp[i+1])) ||
((dedx>tmp[i]) && (dedx<tmp[i+1])))
{
bg[nbg]=(i+1)*0.1 + 0.1*(dedx-tmp[i+1])/(tmp[i]-tmp[i+1]);
nbg++;
if(nbg>2)
{
printf("mTecPIDModule(GetBetaGamma) -> ERROR: Bad nbg %d\n",nbg);
nbg=2;
}
}
}
*nbtgm = nbg;
for(i=0; i<nbg; i++)
{
*(betagamma_ptr+i)=bg[i];
}
return 0;
}
//=====================================================================
int mTecPIDModule::GetMomentum(int itrk, int Debug,
dCglTrackWrapper* dCglTrack,
dCglParticleWrapper* dCglParticle,
float* Momentum, int* dchtrkid) {
// Finds momentum from global tracking
float px,py,pz;
int i,j,found1,found2;
*Momentum = -1.;
*dchtrkid = -1;
found1=0; found2=0;
for(i=0; i<(int)dCglTrack->RowCount(); i++)
{
if(itrk==dCglTrack->get_tectrackid(i) &&
dCglTrack->get_dctracksid(i)>=0)
{
found1=1;
*dchtrkid = dCglTrack->get_dctracksid(i);
for(j=0; j<(int)dCglParticle->RowCount(); j++)
{
if(dCglParticle->get_trackid(j)==dCglTrack->get_id(i))
{
found2=1;
px = dCglParticle->get_pxyz(0,j);
py = dCglParticle->get_pxyz(1,j);
pz = dCglParticle->get_pxyz(2,j);
*Momentum = sqrt(px*px+py*py+pz*pz);
break;
}
} // j
}
} // i
if(*Momentum<0.) return -1;
return 0;
} // end of GetMomentum
//===========================================================================
int mTecPIDModule::GetPerfectMomentum(int itrk, PHCompositeNode *root,
int gtrknum,
float* Momentum, int* Particle, int Debug) {
// Finds Geant momentum for Tec track itrk.
// Closest (in space) Geant track is selected.
int firstplane=0; int lastplane=3;
int nplanes=lastplane-firstplane+1;
PHNodeIterator i(root);
//------ Get Tables ---------------------------------------------
fkinNode_t* fkN = static_cast<fkinNode_t*>(i.findFirst("PHIODataNode", "fkin"));
if (!fkN) {
cerr << "mTecPIDModule -> ERROR: fkin table not found." << endl;
return False;
}
fkinWrapper* fkin = fkN->getData();
//----------------------------------------------------------------
// TecOut
PHTypedNodeIterator<TecOutV1> teciter(root);
TecOutNode_t *TecOutNode = teciter.find("TecOutV1");
TecOutV1* tecout;
if(TecOutNode) {
tecout = (TecOutV1*)TecOutNode->getData();
if(Verbose>0)
cout << "mTecPIDModule: Found TecOut." << endl;
}
else {
if(Verbose>0) cerr << "mTecPIDModule ERROR: Can not find TecOut." << endl;
return False;
}
//-----------------------------------------------------------------
// For every Tec track find closest Geant track
int j,k,jk,trk,gtrack;
float maxdist,xmc[6],ymc[6],chi2,dist[4],x1,x2,y1,y2,Slope = 0.0,Intercept;
j=itrk;
maxdist=99999.;
gtrack=-1;
for(k=0; k<gtrknum; k++) {
for(jk=0; jk<nplanes; jk++) {
xmc[jk]=gtrklist[k].xinglo[jk];
ymc[jk]=gtrklist[k].yinglo[jk];
}
// find distance between Tec track and geant hits
x1 = tecout->getTrackXin(j);
x2 = tecout->getTrackXout(j);
y1 = tecout->getTrackYin(j);
y2 = tecout->getTrackYout(j);
if((x2-x1) != 0.) {
Slope = (y2-y1)/(x2-x1);
Intercept = y1 - Slope*x1;
}
else {
Intercept = 99999.;
}
chi2=0.;
for(jk=0; jk<nplanes; jk++) {
// calculate distance between a point and a line
dist[jk] = (Slope*xmc[jk]-ymc[jk]+Intercept)/sqrt(1.0+Slope*Slope);
chi2 += fabs(dist[jk]);
}
chi2=chi2/((float)(nplanes));
if(chi2<maxdist) {
gtrack=k; maxdist=chi2;
}
} // k loop over geant tracks
*Particle=-1;
*Momentum=-1.;
trk=gtrklist[gtrack].true_track;
for(k=0; k<(int)fkin->RowCount(); k++) {
if(trk==fkin->get_true_track(k)) {
*Particle=fkin->get_idpart(k);
*Momentum=fkin->get_ptot(k);
break;
}
}
if(Debug>0) cout << "trackm: " << j << " " << gtrack << " " << maxdist <<
" " << *Momentum << " " << *Particle << endl;
return 0;
}
//===========================================================================
int mTecPIDModule::CalculateDZ(PHCompositeNode *root, int itrk, int Debug,
int *pc1id, int *pc3id, float *DZ) {
// Calculates DZ = tan(theta) using Z from Pc1 and Pc3
// Association cuts TecPc1Cut and TecPc3 cut are data members
*DZ=-999.;
*pc1id = -1;
*pc3id = -1;
int found=0;
// Get pointers to the tables
PHTypedNodeIterator<dPadClusterWrapper> iPC1N(root);
PHTypedNodeIterator<dPadClusterWrapper> iPC3N(root);
dPadClusterNode_t *PC1N = iPC1N.find("dPc1Cluster");
dPadClusterNode_t *PC3N = iPC3N.find("dPc3Cluster");
if(!PC1N || !PC3N) {
if(Verbose>1)
cerr << "mTecPIDModule::CalculateDZ ERROR: PadCluster tables not found." << endl;
return 1;
}
dPadClusterWrapper* dPc1Cluster = PC1N->getData();
dPadClusterWrapper* dPc3Cluster = PC3N->getData();
// TecOut
PHTypedNodeIterator<TecOutV1> teciter(root);
TecOutNode_t *TecOutNode = teciter.find("TecOutV1");
TecOutV1* tecout;
if(TecOutNode) {
tecout = (TecOutV1*)TecOutNode->getData();
}
else {
if(Verbose>0) cerr << "mTecPIDModule ERROR: Can not find TecOut." << endl;
return -1;
}
// Associate Tec track with Pc1 and Pc3 clusters
float in[3],out[3];
in[0]=tecout->getTrackXin(itrk);
in[1]=tecout->getTrackYin(itrk);
in[2]=-1.;
if(tecout->getTrackSide(itrk)==1) in[2]=1;
out[0]=tecout->getTrackXout(itrk);
out[1]=tecout->getTrackYout(itrk);
out[2]=-1;
if(tecout->getTrackSide(itrk)==1) out[2]=1;
TecTrack* tmp = new TecTrack(in,out);
tmp->project2PC(root);
// Try to find Z from Pc1 and Pc3
if(tmp->getPc1pointer()>-1 && tmp->getPc1distance()<TecPc1Cut &&
tmp->getPc3pointer()>-1 && tmp->getPc3distance()<TecPc3Cut) {
*pc1id = tmp->getPc1pointer();
*pc3id = tmp->getPc3pointer();
float pc1x = dPc1Cluster->get_xyz(0,*pc1id);
float pc1y = dPc1Cluster->get_xyz(1,*pc1id);
float pc1z = dPc1Cluster->get_xyz(2,*pc1id);
float pc3x = dPc3Cluster->get_xyz(0,*pc3id);
float pc3y = dPc3Cluster->get_xyz(1,*pc3id);
float pc3z = dPc3Cluster->get_xyz(2,*pc3id);
float pc1r = sqrt(pc1x*pc1x+pc1y*pc1y);
float pc3r = sqrt(pc3x*pc3x+pc3y*pc3y);
if((pc3r-pc1r)!=0.) {
*DZ = (pc3z-pc1z)/(pc3r-pc1r);
found=1;
}
}
delete tmp;
if(found==1) return 0; else return -1;
}
//===========================================================================
int mTecPIDModule::GetTrackLengthFromPAD(int itrk,
PHCompositeNode *root,
int Debug,
int *pc1id, int *pc3id,
float* TrackLength) {
// Calculates track length using X and Y from Tec and Z from Pc1 and Pc3
*TrackLength = -999.;
*pc1id = -1;
*pc3id = -1;
// Get pointers to the tables
//----------------------------------------------------------------
// TecOut
PHTypedNodeIterator<TecOutV1> teciter(root);
TecOutNode_t *TecOutNode = teciter.find("TecOutV1");
TecOutV1* tecout;
if(TecOutNode) {
tecout = (TecOutV1*)TecOutNode->getData();
}
else {
if(Verbose>0) cerr << "mTecPIDModule ERROR: Can not find TecOut." << endl;
return -1;
}
// TGO
TecGeometryObject *TGO;
ObjectNode_t *TGONode;
PHNodeIterator nodeIter(root);
PHNodeIterator *parNodeIter;
PHCompositeNode *parNode;
parNode =
static_cast < PHCompositeNode * >(nodeIter.findFirst ("PHCompositeNode", "PAR"));
if (!parNode) {
PHMessage ("mTecPIDModule: ", PHError,
"PAR node does not exist.");
return -1;
}
parNodeIter = new PHNodeIterator(parNode);
TGONode =
static_cast < ObjectNode_t * >(parNodeIter->findFirst ("PHIODataNode", "TecGeometry"));
if (!TGONode) {
PHMessage ("mTecPIDModule: ", PHError,
"TecGeometryObject not found.\n");
delete parNodeIter;
return -1;
}
else {
TGO = static_cast < TecGeometryObject * >(TGONode->getData ());
}
delete parNodeIter;
//---------------------------------------------------------------------------
float beg1,beg2;
mTecPIDModule::GetBeg12(TGO, &beg1, &beg2);
// Use dx and dy from Tec track
float dx = tecout->getTrackXout(itrk)-tecout->getTrackXin(itrk);
float dy = tecout->getTrackYout(itrk)-tecout->getTrackYin(itrk);
float dz;
float DZ=-999.;
int found=0;
// Try to find Z from Pc1 and Pc3
int statdz = CalculateDZ(root, itrk, Debug, pc1id, pc3id, &DZ);
if(statdz==0 && DZ!=-999.) {
dz = fabs(beg2-beg1)*DZ;
*TrackLength = sqrt(dx*dx+dy*dy+dz*dz);
found=1;
}
if(found==1) return 0; else return -1;
}
//========================================================================
int mTecPIDModule::GetTrackLength(int itrk,
TecOutV1* tecout,
dCglTrackWrapper* dCglTrack,
dDchTracksWrapper* dDchTracks,
float& TrackLength) {
// Calculates track length using X and Y from Tec and Z from Dch
float dx = tecout->getTrackXout(itrk)-tecout->getTrackXin(itrk);
float dy = tecout->getTrackYout(itrk)-tecout->getTrackYin(itrk);
float trl = sqrt(dx*dx+dy*dy);
float dz = 0.0;
// find corresponding Dch track (if any)
int found=0;
for(int i=0; i<(int)dCglTrack->RowCount(); i++) {
if(itrk==dCglTrack->get_tectrackid(i) && dCglTrack->get_dctracksid(i)>=0) {
int dchtrkid = dCglTrack->get_dctracksid(i);
float dxtmp = dDchTracks->get_direction(0,dchtrkid);
float dytmp = dDchTracks->get_direction(1,dchtrkid);
float dztmp = dDchTracks->get_direction(2,dchtrkid);
if((fabs(dxtmp)+fabs(dytmp))!=0) {
dz = dztmp/sqrt(dxtmp*dxtmp+dytmp*dytmp);
found=1;
break;
}
}
}
if(found) {
TrackLength = trl*sqrt(1.0+dz*dz);
return 1;
}
else {
TrackLength = -trl;
return -1;
}
}
//========================================================================
int mTecPIDModule::GetTrackLength(int itrk, int Debug,
TecOutV1* tecout,
TecGeometryObject* TGO,
dCglTrackWrapper* dCglTrack,
dDchTracksWrapper* dDchTracks,
int *dchtrkid, float* TrackLength) {
// Calculates track length using X and Y from Tec and Z from Dch
*TrackLength = 0.;
*dchtrkid = -1;
float dx = tecout->getTrackXout(itrk)-tecout->getTrackXin(itrk);
float dy = tecout->getTrackYout(itrk)-tecout->getTrackYin(itrk);
float dz = 0.0;
int found=0;
float beg1,beg2;
mTecPIDModule::GetBeg12(TGO, &beg1, &beg2);
// find corresponding Dch track (if any)
for(int i=0; i<(int)dCglTrack->RowCount(); i++) {
if(itrk==dCglTrack->get_tectrackid(i) && dCglTrack->get_dctracksid(i)>=0) {
*dchtrkid = dCglTrack->get_dctracksid(i);
float dxtmp = dDchTracks->get_direction(0,*dchtrkid);
float dytmp = dDchTracks->get_direction(1,*dchtrkid);
float dztmp = dDchTracks->get_direction(2,*dchtrkid);
if((fabs(dxtmp)+fabs(dytmp))!=0) {
dz = dztmp/sqrt(dxtmp*dxtmp+dytmp*dytmp);
found=1;
break;
}
}
}
if(found==1) {
dz = fabs(beg2-beg1)*dz;
*TrackLength = sqrt(dx*dx+dy*dy+dz*dz);
if(Debug > 1) cout << "mTecPID -> itrk,TrkLen: " << itrk << " "
<< *TrackLength << endl;
return 0;
}
else {
if(Debug > 1) cout << "mTecPID -> itrk,TrkLen: " << itrk << " "
<< *TrackLength << endl;
return 1;
}
}
//===========================================================================
int mTecPIDModule::GetPerfectTrackLength(int itrk, PHCompositeNode *root,
int gtrknum,
float* TrackLength, int Debug) {
int firstplane=0; int lastplane=3;
int nplanes=lastplane-firstplane+1;
PHNodeIterator i(root);
//------ Get Tables ---------------------------------------------
tecghitNode_t* tgN = static_cast<tecghitNode_t*>(i.findFirst("PHIODataNode", "tecghit"));
if (!tgN) {
cerr << "mTecPIDModule -> ERROR: tecghit table not found." << endl;
return False;
}
tecghitWrapper* tecghit = tgN->getData();
//----------------------------------------------------------------
// TecOut
PHTypedNodeIterator<TecOutV1> teciter(root);
TecOutNode_t *TecOutNode = teciter.find("TecOutV1");
TecOutV1* tecout;
if(TecOutNode) {
tecout = (TecOutV1*)TecOutNode->getData();
}
else {
if(Verbose>0) cerr << "mTecPIDModule ERROR: Can not find TecOut." << endl;
return -1;
}
// TGO
TecGeometryObject *TGO;
ObjectNode_t *TGONode;
PHNodeIterator nodeIter(root);
PHNodeIterator *parNodeIter;
PHCompositeNode *parNode;
parNode =
static_cast < PHCompositeNode * >(nodeIter.findFirst ("PHCompositeNode", "PAR"));
if (!parNode) {
PHMessage ("mTecPIDModule: ", PHError,
"PAR node does not exist.");
return -1;
}
parNodeIter = new PHNodeIterator(parNode);
TGONode =
static_cast < ObjectNode_t * >(parNodeIter->findFirst ("PHIODataNode", "TecGeometry"));
if (!TGONode) {
PHMessage ("mTecPIDModule: ", PHError,
"TecGeometryObject not found.\n");
delete parNodeIter;
return -1;
}
else {
TGO = static_cast < TecGeometryObject * >(TGONode->getData ());
}
delete parNodeIter;
//---------------------------------------------------------------------------
// For every Tec track find closest Geant track
*TrackLength=-999.;
int j,k,jk,trk,gtrack;
float maxdist,xmc[6],ymc[6];
j=itrk;
maxdist=99999.;
gtrack=-1;
for(k=0; k<gtrknum; k++) {
for(jk=0; jk<nplanes; jk++) {
xmc[jk]=gtrklist[k].xinglo[jk];
ymc[jk]=gtrklist[k].yinglo[jk];
}
// find distance between Tec track and geant hits
float chi2,dist[4],x1,x2,y1,y2,Slope = 0.0,Intercept;
x1 = tecout->getTrackXin(j);
x2 = tecout->getTrackXout(j);
y1 = tecout->getTrackYin(j);
y2 = tecout->getTrackYout(j);
if((x2-x1) != 0.)
{
Slope = (y2-y1)/(x2-x1);
Intercept = y1 - Slope*x1;
}
else {
Intercept = 99999.;
}
chi2=0.;
for(jk=0; jk<nplanes; jk++)
{
// calculate distance between a point and a line
dist[jk] = (Slope*xmc[jk]-ymc[jk]+Intercept)/sqrt(1.0+Slope*Slope);
chi2 += fabs(dist[jk]);
}
chi2=chi2/((float)(nplanes));
if(chi2<maxdist) {
gtrack=k; maxdist=chi2;
}
} // k loop over geant tracks
// find track length
float dx = tecout->getTrackXout(itrk)-tecout->getTrackXin(itrk);
float dy = tecout->getTrackYout(itrk)-tecout->getTrackYin(itrk);
float beg1,beg2;
mTecPIDModule::GetBeg12(TGO, &beg1, &beg2);
float x1,x2,y1,y2,z1,z2,zdirection = 0.0;
int found=0;
trk=gtrklist[gtrack].true_track;
for(k=0; k<(int)tecghit->RowCount(); k++) {
if(trk==tecghit->get_mctrack(k)) {
// local tec coordinate system is different from global.
x1=tecghit->get_xyzinloc(0,k);
x2=tecghit->get_xyzoutloc(0,k);
y1=tecghit->get_xyzinloc(2,k);
y2=tecghit->get_xyzoutloc(2,k);
z1=tecghit->get_xyzinloc(1,k);
z2=tecghit->get_xyzoutloc(1,k);
zdirection = fabs(z2-z1)/fabs(y2-y1);
found=1;
if(Debug>1) cout << "xyz: " << x1 << " " << x2 << " " << y1 << " " << y2
<< " " << z1 << " " << z2 << " " << endl;
break;
}
}
if(found==1) {
float dz = fabs(beg2-beg1)*zdirection;
*TrackLength = sqrt(dx*dx+dy*dy+dz*dz);
}
if(Debug>0) cout << "trackl: " << j << " " << gtrack << " " << maxdist
<< " " << *TrackLength << endl;
return 0;
}
//===========================================================================
int mTecPIDModule::GuessMomentum(int itrk, PHCompositeNode *root,
float* Momentum, float* Charge) {
// Calculates transverse momentum using Tec Alpha angle
//----------------------------------------------------------------
// TecOut
PHTypedNodeIterator<TecOutV1> teciter(root);
TecOutNode_t *TecOutNode = teciter.find("TecOutV1");
TecOutV1* tecout;
if(TecOutNode) {
tecout = (TecOutV1*)TecOutNode->getData();
}
else {
if(Verbose>0) cerr << "mTecPIDModule ERROR: Can not find TecOut." << endl;
return -1;
}
*Momentum = tecout->getTrackPt(itrk);
*Charge = tecout->getTrackCharge(itrk);
return 0;
}
//===========================================================================
mTecPIDModule::~mTecPIDModule(){}
//===========================================================================
//===========================================================================
//===========================================================================
mTecPIDModule::mTecPIDModule(){
Verbose=0;
method=1;
Truncation=0.6;
mip=5.35;
el_gaindep[0]=1.71;
el_gaindep[1]=-0.0212;
el_gaindep[2]=0.;
mip_gaindep[0]=-0.00484;
mip_gaindep[1]=0.1014;
mip_gaindep[2]=0.;
el_trnkdep[0]=1.70;
el_trnkdep[1]=-0.348;
el_trnkdep[2]=0.;
mip_trnkdep[0]=0.914;
mip_trnkdep[1]=-2.83;
mip_trnkdep[2]=4.95;
Sigma[0]=0.20;
Sigma[1]=0.15;
Sigma[2]=0.15;
Sigma[3]=0.15;
numbinmin=100;
Calibration[0]=1.41;
Calibration[1]=10000.0;
UseOnlyMerged=0;
Write2File=0;
Write2Ntuple=0;
RejectOverlaps=0; // overlaps are always rejected
PerfectPID=0;
UseTec=0;
TecPc1Cut = 2.5;
TecPc3Cut = 4.0;
hitcut=3;
}
//===========================================================================
PHBoolean mTecPIDModule::analysesim(TChain* chain1) {
using namespace TecUtilities;
TFile* OutputNtupleFile = new TFile("test.root","RECREATE");
TNtuple* ntpcheck = new TNtuple("ntpcheck"," ","evt:nhits:ntrunk:trnksum:dphiproj:dphipc3:mom:toftof:tofpl:bbct0:charge:tecalpha:dchalpha:cent:tecmult:dchmult");
float ntp[16];
vector<float> buffer;
multimap<int, float, less<int> > HITMAP;
rng = gsl_rng_alloc(gsl_rng_mt19937);
// PadCluster* pc3 = 0;
dPadClusterWrapper* pc3 = 0;
CglTrack* cgl = 0;
// dCglTrackWrapper* cgl = 0;
dTofReconstructedWrapper* tof = 0;
DchTrack* dch = 0;
TecOut* tec = 0;
PHTrackOut* proj = 0;
chain1->SetBranchAddress("DST/dPc3Cluster", &pc3);
chain1->SetBranchAddress("DST/DchTrack", &dch);
chain1->SetBranchAddress("DST/CglTrack", &cgl);
chain1->SetBranchAddress("DST/dTofReconstructed", &tof);
chain1->SetBranchAddress("DST/TecOut", &tec);
chain1->SetBranchAddress("DST/PHTrackOut", &proj);
cout << "Total events: " << chain1->GetEntries() << endl;
// chain1->Print();
// loop over events
for(int i=1; i<chain1->GetEntries(); i++) {
if(i%100==0) cout << "event # " << i << endl;
chain1->GetEvent(i);
// make hit list
for(int j=0; j<tec->getNHits(); j++) {
int trkid=tec->getHitTrackID(j);
if(trkid>-1) {
// int jindex = tec->getHitIndex(j);
int myadc = tec->getHitADC(j);
float mycharge = Ampl2Charge(myadc);
float gainA = 1.;
mycharge = mycharge*gainA;
if(myadc>2) HITMAP.insert(pair<int, float>(trkid,mycharge));
}
}
// cout << i << " " << tec->get_TecNTrack() << " " << tec->getNHits() << " " << cgl->get_CglNTrack()
// << " " << tof->RowCount() << " " << pc3->RowCount() << " " << dch->get_DchNTrack() << endl;
for(int itrk=0; itrk<(int)tec->get_TecNTrack(); itrk++) {
float dchalpha = -999.;
float tecalpha = tec->getTrackAlpha(itrk);
float tecphi = tec->getTrackPhi(itrk);
float tecnhits = (float)tec->getTrackNhits(itrk);
// int tecsect = tec->getTrackSector(itrk);
float pc3phi = 9999.;
float projphi = 9999.;
float dphiproj = 9999.;
float dphipc3 = 9999.;
float mom = 0.;
float toftof = 9999.;
float tofpl = 999.;
float charge = 1.;
// if(tecsect==1) {
// Check if this Tec track was associated with a global track
// and calculate phi difference for Dch track projection and Pc3 cluster
for(int icgl=0; icgl<(int)cgl->get_CglNTrack(); icgl++) {
if(cgl->get_tectrackid(icgl)==itrk) {
float xp = proj->get_projectionTec(icgl,0);
float yp = proj->get_projectionTec(icgl,1);
if(xp!=0) projphi = 3.1415927+atan(yp/xp);
if(projphi<0) projphi=2*3.1415927+projphi;
dphiproj = tecphi - projphi;
int pc3ptr = cgl->get_pc3clusid(icgl);
if(pc3ptr>-1 && pc3ptr<(int)pc3->RowCount()) {
float x = pc3->get_xyz(0,pc3ptr);
float y = pc3->get_xyz(1,pc3ptr);
if(x!=0) pc3phi = 3.1415927+atan(y/x);
if(pc3phi<0) pc3phi=2*3.1415927+pc3phi;
dphipc3 = tecphi-pc3phi;
}
if(dch->get_alpha(icgl)>0) charge = -1.;
mom = dch->get_momentum(icgl);
dchalpha = dch->get_alpha(icgl);
int tofptr = cgl->get_tofrecid(icgl);
if(tofptr>-1) {
toftof = tof->get_tof(tofptr);
}
tofpl = proj->get_tofPathLength(icgl);
break;
}
}
if(toftof<999.) {
// cout << " good track: " << mom << " " << toftof << " " << tofpl << " " << tecphi << " " << dphiproj <<
// " " << dphipc3 << endl;
// Calculate truncated dE/dX
// Fill buffer from hitmap
for(multimap<int, float, less<int> >::iterator it = HITMAP.begin(); it != HITMAP.end(); ++it)
{
if((*it).first==itrk) {
buffer.push_back((*it).second);
}
}
// Re-order buffer
sort(buffer.begin(), buffer.end());
float truncation = 0.6;
int ntrunk = ((int)(buffer.size()*truncation));
float trnksum = 0.;
for(int k=0; k<ntrunk; k++) trnksum += buffer[k];
ntp[0]=(float)i;
ntp[1]=(float)tecnhits;
ntp[2]=ntrunk;
ntp[3]=trnksum;
ntp[4]=dphiproj;
ntp[5]=dphipc3;
ntp[6]=mom;
ntp[7]=toftof;
ntp[8]=tofpl;
ntp[9]=0.;
ntp[10]=charge;
ntp[11]=tecalpha;
ntp[12]=dchalpha;
ntp[13]=0.;
ntp[14]=(float)tec->get_TecNTrack();
ntp[15]=(float)cgl->get_CglNTrack();
ntpcheck->Fill(ntp);
} // good tof
// } // sector 1
buffer.clear();
} // end loop over tec tracks
HITMAP.clear();
} // end loop over events
OutputNtupleFile->Write();
OutputNtupleFile->Close();
return True;
}
//===========================================================================
PHBoolean mTecPIDModule::myglobaltracking(TChain* chain1, TChain* chain2, TecCalibrationObject* TCO) {
using namespace PhUtilities;
using namespace TecUtilities;
//chain1->Print();
TFile* OutputNtupleFile = new TFile("test.root","RECREATE");
TNtuple* ntpcheck = new TNtuple("ntpcheck"," ","evt:nhits:ntrunk:trnksum:dphiproj:dphipc3:mom:toftof:tofpl:bbct0:charge:tecalpha:dchalpha:cent:tecmult:dchmult:dphiprojn:dalphan:distn");
float ntp[19];
float tecalpha0 = 14.1 / 1000.0; // parameters for calculating sigma Alpha
float tecalpha1 = 101.4 / 1000.0;
// float tecalpha2 = 0.0 / 1000.0;
float tecphi0 = 0.754 / 1000.0; // parameters for calculating sigma Phi
float tecphi1 = 15.4 / 1000.0;
// float tecphi2 = 0.0 / 1000.0;
double a0 = tecalpha0 / 2.;
double a1 = tecalpha1 / 2.;
double b0 = tecphi0;
double b1 = tecphi1;
vector<float> buffer;
multimap<int, float, less<int> > HITMAP;
rng = gsl_rng_alloc(gsl_rng_mt19937);
EventHeader* evt = 0;
RunHeader* run = 0;
// PadCluster* pc3 = 0;
dPadClusterWrapper* pc3 = 0;
// CglTrack* cgl = 0;
dCglTrackWrapper* cgl = 0;
dTofReconstructedWrapper* tof = 0;
DchTrack* dch = 0;
TecOut* tec = 0;
PHTrackOut* proj = 0;
BbcOut* bbc = 0;
ZdcOut* zdc = 0;
chain2->SetBranchAddress("RUN/RunHeader", &run);
chain1->SetBranchAddress("DST/EventHeader", &evt);
chain1->SetBranchAddress("DST/dPc3Cluster", &pc3);
chain1->SetBranchAddress("DST/DchTrack", &dch);
chain1->SetBranchAddress("DST/dCglTrack", &cgl);
chain1->SetBranchAddress("DST/dTofReconstructed", &tof);
chain1->SetBranchAddress("DST/TecOut", &tec);
chain1->SetBranchAddress("DST/PHTrackOut", &proj);
chain1->SetBranchAddress("DST/BbcOut", &bbc);
chain1->SetBranchAddress("DST/ZdcOut", &zdc);
// run info
chain2->GetEvent(1);
cout << "Total events: " << chain1->GetEntries() << endl;
// loop over events
// for(int i=1; i<chain1->GetEntries(); i++) {
for(int i=1; i<200; i++) {
if(i%1==0) cout << "---------- event # " << i << endl;
chain1->GetEvent(i);
// make hit list
for(int j=0; j<tec->getNHits(); j++) {
int trkid=tec->getHitTrackID(j);
if(trkid>-1) {
int jindex = tec->getHitIndex(j);
int myadc = tec->getHitADC(j);
float mycharge = Ampl2Charge(myadc);
float gainA = TCO->getAbsoluteGain(jindex);
mycharge = mycharge*gainA;
HITMAP.insert(pair<int, float>(trkid,mycharge));
}
}
cout << " tec: " << tec->get_TecNTrack() << endl;
float bbct0 = bbc->get_TimeZero();
float bbc1 = bbc->get_ChargeSum(Bbc::South);
float bbc2 = bbc->get_ChargeSum(Bbc::North);
float zdc1 = zdc->get_Energy(Zdc::South);
float zdc2 = zdc->get_Energy(Zdc::North);
int cent = getCentralityByClock(bbc1, bbc2, zdc1, zdc2);
for(int icgl=0; icgl<(int)cgl->RowCount(); icgl++) {
if(dch->get_arm(icgl)==0) {
float pc3phi = 9999.;
float projphi = 9999.;
float projphipc3 = 9999.;
float dphiproj = 9999.;
float dalpha = 9999.;
float dphipc3 = 9999.;
float mom = 0.;
float toftof = 9999.;
float tofpl = 999.;
float charge = 1.;
float xp = proj->get_projectionTec(icgl,0);
float yp = proj->get_projectionTec(icgl,1);
if(xp!=0) projphi = 3.1415927+atan(yp/xp);
if(projphi<0) projphi=2*3.1415927+projphi;
float xppc3 = proj->get_projectionPc3(icgl,0);
float yppc3 = proj->get_projectionPc3(icgl,1);
if(xppc3!=0) projphipc3 = 3.1415927+atan(yppc3/xppc3);
if(projphi<0) projphi=2*3.1415927+projphi;
float dchalpha = dch->get_alpha(icgl);
float dchphi = dch->get_phi(icgl);
float dchphi0 = dch->get_phi0(icgl);
int pc3ptr = cgl->get_pc3clusid(icgl);
if(pc3ptr>-1 && pc3ptr<(int)pc3->RowCount()) {
float x = pc3->get_xyz(0,pc3ptr);
float y = pc3->get_xyz(1,pc3ptr);
if(x!=0) pc3phi = 3.1415927+atan(y/x);
if(pc3phi<0) pc3phi=2*3.1415927+pc3phi;
dphipc3 = projphi-pc3phi;
}
if(dch->get_alpha(icgl)>0) charge = -1.;
mom = dch->get_momentum(icgl);
int tofptr = cgl->get_tofrecid(icgl);
if(tofptr>-1) {
toftof = tof->get_tof(tofptr);
}
tofpl = proj->get_tofPathLength(icgl);
cout << icgl << ": " << projphi << " " << projphipc3 << " "
<< dchphi << " " << dchphi0 << " " << pc3phi << " "
<< mom << " " << dchalpha << endl;
int tecptr = -1;
float dist=9999.;
float mindist=9999.;
float distfound=9999.;
float dphiprojfound=9999.;
float dalphafound=9999.;
float tecalpha=999.;
float tecphi=999.;
float tecnhits=0.;
int tecsect = 99;
double sigmaalpha = a0 + a1 * fabs(dchalpha);
double sigmaphi = b0 + b1 * fabs(dchalpha);
// Loop over tec tracks
for(int itrk=0; itrk<(int)tec->get_TecNTrack(); itrk++) {
tecalpha = tec->getTrackAlpha(itrk);
tecphi = tec->getTrackPhi(itrk);
tecnhits = (float)tec->getTrackNhits(itrk);
tecsect = tec->getTrackSector(itrk);
// Find closest TEC track
dphiproj = (tecphi-projphi)/sigmaphi;
dalpha = (tecalpha*2.-dchalpha)/sigmaalpha;
dist = sqrt(dphiproj*dphiproj+dalpha*dalpha);
if(dist<mindist) {
mindist=dist;
tecptr=itrk;
dphiprojfound = dphiproj;
dalphafound = dalpha;
distfound = dist;
}
} // end loop over tec tracks
// Calculate truncated dE/dX for closest Tec track
if(tecptr>-1) {
cout << " tec: " << tec->getTrackPhi(tecptr) << " "
<< tec->getTrackAlpha(tecptr) << " " << dchalpha << " "
<< sigmaalpha << " " << sigmaphi << " " << dphiprojfound << " "
<< dalphafound << " " << distfound << endl;
// Fill buffer from hitmap
for(multimap<int, float, less<int> >::iterator it = HITMAP.begin(); it != HITMAP.end(); ++it)
{
if((*it).first==tecptr) {
buffer.push_back((*it).second);
}
}
// Re-order buffer
sort(buffer.begin(), buffer.end());
float truncation = 0.6;
int ntrunk = ((int)(buffer.size()*truncation));
float trnksum = 0.;
for(int k=0; k<ntrunk; k++) trnksum += buffer[k];
ntp[0]=(float)i;
ntp[1]=(float)tecnhits;
ntp[2]=ntrunk;
ntp[3]=trnksum;
ntp[4]=dphiproj;
ntp[5]=dphipc3;
ntp[6]=mom;
ntp[7]=toftof;
ntp[8]=tofpl;
ntp[9]=bbct0;
ntp[10]=charge;
ntp[11]=tecalpha;
ntp[12]=dchalpha;
ntp[13]=(float)cent;
ntp[14]=(float)tec->get_TecNTrack();
ntp[15]=(float)cgl->RowCount();
ntp[16]=dphiprojfound;
ntp[17]=dalphafound;
ntp[18]=distfound;
ntpcheck->Fill(ntp);
buffer.clear();
} // if tecptr>-1
} // east arm
} // end loop over glogal tracks
HITMAP.clear();
} // end loop over events
OutputNtupleFile->Write();
OutputNtupleFile->Close();
return True;
}
//===========================================================================
void mTecPIDModule::openfile(char* fname) {
OutputNtupleFile = new TFile(fname,"RECREATE");
cout << "output file " << fname << " opened." << endl;
ntpcheck = new TNtuple("ntpcheck"," ","evt:nhits:ntrunk:trnksum:dphiproj:dphipc3:mom:toftof:tofpl:bbct0:charge:tecalpha:dchalpha:cent:tecmult:dchmult:trnksum05:trnksum07:trnksum10:ntrunk05:ntrunk07:ntrunk10:trklen:el:npmt0:ntrunk08:ntrunk09:trnksum08:trnksum09:pc3z:projz:zcorr:nhits0:nhits1:nhits2:nhits3:nhits4:nhits5:index");
}
//===========================================================================
void mTecPIDModule::closefile() {
cout << "writing out ntuple..." << endl;
OutputNtupleFile->Write();
cout << "closing file..." << endl;
OutputNtupleFile->Close();
cout << "file closed." << endl;
}
float get_phiv(float px1, float py1, float pz1, float px2, float py2, float pz2) {
float phiv=999.;
float ux = px1+px2; // normalized photon direction
float uy = py1+py2;
float uz = pz1+pz2;
float uu = sqrt(ux*ux+uy*uy+uz*uz);
ux=ux/uu;
uy=uy/uu;
uz=uz/uu;
float vx = py1*pz2-pz1*py2; // normalized perpendicular to decay plane
float vy = -px1*pz2+pz1*px2;
float vz = px1*py2-py1*px2;
float vv = sqrt(vx*vx+vy*vy+vz*vz);
vx=vx/vv;
vy=vy/vv;
vz=vz/vv;
float wx = uy*vz-uz*vy; // perpendicular to both u and v
float wy = -ux*vz+uz*vx;
float wz = ux*vy-uy*vx;
float uax = uy/sqrt(uy*uy+ux*ux);
float uay = -ux/sqrt(uy*uy+ux*ux);
float uaz = 0.;
phiv = acos(wx*uax+wy*uay+wz*uaz); // angle between mag. field direction and
// perpendicular to decay plane
return phiv;
}
float get_invmass(float px1, float py1, float pz1, float px2, float py2, float pz2) {
float invmass=0.;
float e1 = sqrt(px1*px1+py1*py1+pz1*pz1);
float e2 = sqrt(px2*px2+py2*py2+pz2*pz2);
float dotproduct = px1*px2+py1*py2+pz1*pz2;
float mass2 = 2.*(e1*e2-dotproduct);
if(mass2>0.) {
invmass=sqrt(mass2);
}
return invmass;
}
//===========================================================================
PHBoolean mTecPIDModule::analysemerged(DstContent* dst, TecCalibrationObject* TCO) {
// BbcOut *bbc = static_cast<BbcOut*> (dst->getClass("BbcOut"));
TecOut *tec = static_cast<TecOut*> (dst->getClass("TecOut"));
DchTrack *dch = static_cast<DchTrack*> (dst->getClass("DchTrack"));
CglTrack* cgl = static_cast<CglTrack*> (dst->getClass("CglTrack"));
TofOut* tof = static_cast<TofOut*> (dst->getClass("TofOut"));
PadCluster* pc3 = static_cast<PadCluster*> (dst->getClass("Pc3Cluster"));
// PHTrackOut* proj = static_cast<PHTrackOut*> (dst->getClass("PHTrackOut"));
cout << "event # " << evtnum << endl;
cout << "$$$$$$$ tec: " << tec->get_TecNTrack() << endl;
cout << "$$$$$$$ dch: " << dch->get_DchNTrack() << endl;
cout << "$$$$$$$ cgl: " << cgl->get_CglNTrack() << endl;
cout << "$$$$$$$ pad: " << pc3->get_PadNCluster() << endl;
cout << "$$$$$$$ tof: " << tof->get_TofNHit() << endl;
evtnum++;
return True;
}
//===========================================================================
PHBoolean mTecPIDModule::analysesim(DstContent* dst, TecCalibrationObject* TCO) {
using namespace PhUtilities;
using namespace TecUtilities;
vector<float> buffer;
multimap<int, float, less<int> > HITMAP;
float ntp[39];
// TriggerHelper triggerhelper(dst);
// L2DecisionHelper lvl2Helper(dst);
// TrigLvl1 *trig = static_cast<TrigLvl1*> (dst->getClass("TrigLvl1"));
// TrigRunLvl1* trigrun1 = static_cast<TrigRunLvl1*> (dst->getClass("TrigRunLvl1"));
// TrigRunLvl2* trigrun2 = static_cast<TrigRunLvl2*> (dst->getClass("TrigRunLvl2"));
// Lvl2DecisionOut* lvl2 = static_cast<Lvl2DecisionOut*> (dst->getClass("L2Decision"));
// RunHeader *run = static_cast<RunHeader*> (dst->getClass("RunHeader"));
// EventHeader *evt = static_cast<EventHeader*> (dst->getClass("EventHeader"));
BbcOut *bbc = static_cast<BbcOut*> (dst->getClass("BbcOut"));
// VtxOut *vtx = static_cast<VtxOut*> (dst->getClass("VtxOut"));
// ZdcOut *zdc = static_cast<ZdcOut*> (dst->getClass("ZdcOut"));
TecOut *tec = static_cast<TecOut*> (dst->getClass("TecOut"));
DchTrack *dch = static_cast<DchTrack*> (dst->getClass("DchTrack"));
CglTrack* cgl = static_cast<CglTrack*> (dst->getClass("CglTrack"));
// dCglTrackWrapper* cgl = static_cast<dCglTrackWrapper*> (dst->getClass("dCglTrack"));
// CrkRing* crk = static_cast<CrkRing*> (dst->getClass("CrkRing"));
// dTofReconstructedWrapper* tof = static_cast<dTofReconstructedWrapper*> (dst->getClass("dTofReconstructed"));
TofOut* tof = static_cast<TofOut*> (dst->getClass("TofOut"));
// PadCluster* pc1 = static_cast<PadCluster*> (dst->getClass("Pc1Cluster"));
// PadCluster* pc2 = static_cast<PadCluster*> (dst->getClass("Pc2Cluster"));
// dPadClusterWrapper* pc3 = static_cast<dPadClusterWrapper*> (dst->getClass("dPc3Cluster"));
PadCluster* pc3 = static_cast<PadCluster*> (dst->getClass("Pc3Cluster"));
// EmcClusterLocalExt* emc =
// static_cast<EmcClusterLocalExt*> (dst->getClass("EmcClusterLocalExt"));
PHTrackOut* proj = static_cast<PHTrackOut*> (dst->getClass("PHTrackOut"));
TecGeometryObject* TGO = 0;
static float beg1=0.;
static float beg2=0.;
if (!init_done) {
cout << "Initializing random number generator..." << endl;
rng = gsl_rng_alloc(gsl_rng_mt19937);
TGO = new TecGeometryObject();
TGO->FetchFromFile();
mTecPIDModule::GetBeg12(TGO, &beg1, &beg2);
delete TGO;
init_done = 1;
cout << "beg12: " << beg1 << " " << beg2 << endl;
}
if(evtnum%1==0) cout << "event # " << evtnum << endl;
// make hit list
for(int j=0; j<tec->getNHits(); j++) {
int trkid=tec->getHitTrackID(j);
if(trkid>-1) {
int jindex = tec->getHitIndex(j);
int myadc = tec->getHitADC(j);
float mycharge = Ampl2Charge(myadc);
float gainA = TCO->getAbsoluteGain(jindex);
mycharge = mycharge*gainA;
HITMAP.insert(pair<int, float>(trkid,mycharge));
}
}
// try to find conversion pairs
vector<int> pairs;
/*
for(int i=0; i<(int)dch->get_DchNTrack()-1; i++) {
for(int j=i+1; j<(int)dch->get_DchNTrack(); j++) {
if(dch->get_arm(i)==0 && dch->get_arm(j)==0) {
if((dch->get_alpha(i)*dch->get_alpha(j)<0)) {
float px1 = dch->get_momentum(i)*sin(dch->get_theta0(i))*cos(dch->get_phi0(i));
float py1 = dch->get_momentum(i)*sin(dch->get_theta0(i))*sin(dch->get_phi0(i));
float pz1 = dch->get_momentum(i)*cos(dch->get_theta0(i));
float px2 = dch->get_momentum(j)*sin(dch->get_theta0(j))*cos(dch->get_phi0(j));
float py2 = dch->get_momentum(j)*sin(dch->get_theta0(j))*sin(dch->get_phi0(j));
float pz2 = dch->get_momentum(j)*cos(dch->get_theta0(j));
float phiv = get_phiv(px1,py1,pz1,px2,py2,pz2);
float invmass = get_invmass(px1,py1,pz1,px2,py2,pz2);
if(fabs(phiv)<0.1 && invmass>0.065 && invmass<0.090) {
pairs.push_back(i); pairs.push_back(j);
}
}
}
}
}
*/
float bbct0 = bbc->get_TimeZero();
// float bbc1 = bbc->get_ChargeSum(Bbc::South);
// float bbc2 = bbc->get_ChargeSum(Bbc::North);
// float zdc1 = zdc->get_Energy(Zdc::South);
// float zdc2 = zdc->get_Energy(Zdc::North);
// int cent = getCentralityByClock(bbc1, bbc2, zdc1, zdc2);
int cent = 0;
for(int itrk=0; itrk<(int)tec->get_TecNTrack(); itrk++) {
float dchalpha = -999.;
float tecalpha = tec->getTrackAlpha(itrk);
float tecphi = tec->getTrackPhi(itrk);
float tecnhits = (float)tec->getTrackNhits(itrk);
float pc3phi = 9999.;
float projphi = 9999.;
float dphiproj = 9999.;
float dphipc3 = 9999.;
float mom = 0.;
float toftof = 9999.;
float tofpl = 999.;
float charge = 1.;
// float trklen = tec->getTrackLength(itrk);
float mytrklen = -999.;
float theta0 = 0.;
int el = 0;
int npmt0 = 0;
int tecindex = tec->getTrackIndex(itrk);
// calculate track length
float dx = tec->getTrackXout(itrk)-tec->getTrackXin(itrk);
float dy = tec->getTrackYout(itrk)-tec->getTrackYin(itrk);
float dz = 0.0;
// Check if this Tec track was associated with a global track
// and calculate phi difference for Dch track projection and Pc3 cluster
for(int icgl=0; icgl<(int)cgl->get_CglNTrack(); icgl++) {
if(cgl->get_tectrackid(icgl)==itrk) {
float xp = proj->get_projectionTec(icgl,0);
float yp = proj->get_projectionTec(icgl,1);
if(xp!=0) projphi = 3.1415927+atan(yp/xp);
if(projphi<0) projphi=2*3.1415927+projphi;
dphiproj = tecphi - projphi;
int pc3ptr = cgl->get_pc3clusid(icgl);
if(pc3ptr>-1 && pc3ptr<(int)pc3->get_PadNCluster()) {
float x = pc3->get_xyz(0,pc3ptr);
float y = pc3->get_xyz(1,pc3ptr);
if(x!=0) pc3phi = 3.1415927+atan(y/x);
if(pc3phi<0) pc3phi=2*3.1415927+pc3phi;
dphipc3 = tecphi-pc3phi;
}
if(dch->get_alpha(icgl)>0) charge = -1.;
mom = dch->get_momentum(icgl);
dchalpha = dch->get_alpha(icgl);
theta0 = dch->get_theta0(icgl);
dz = fabs(beg1-beg2)*tan(3.1415927/2.-theta0);
int tofptr = cgl->get_tofrecid(icgl);
if(tofptr>-1) {
toftof = tof->get_tof(tofptr);
tofpl = proj->get_tofPathLength(icgl);
}
// int crkptr = cgl->get_richringid(icgl);
// if(crkptr>-1) {
// npmt0 = crk->get_npmt0(crkptr);
// }
// // check if this is an electron
// for(int i=0; i<(int)pairs.size(); i++) {
// if(icgl==pairs[i] && npmt0>0) el = 1;
// }
break;
}
}
mytrklen = sqrt(dx*dx+dy*dy+dz*dz);
// if(toftof<999.) {
// Calculate truncated dE/dX
// Fill buffer from hitmap
for(multimap<int, float, less<int> >::iterator it = HITMAP.begin(); it != HITMAP.end(); ++it)
{
if((*it).first==itrk) {
buffer.push_back((*it).second);
}
}
// Re-order buffer
sort(buffer.begin(), buffer.end());
int ntrunk05 = (int)(buffer.size()*0.5);
int ntrunk06 = (int)(buffer.size()*0.6);
int ntrunk07 = (int)(buffer.size()*0.7);
int ntrunk08 = (int)(buffer.size()*0.8);
int ntrunk09 = (int)(buffer.size()*0.9);
int ntrunk10 = (int)buffer.size();
float trnksum05 = 0.;
float trnksum06 = 0.;
float trnksum07 = 0.;
float trnksum08 = 0.;
float trnksum09 = 0.;
float trnksum10 = 0.;
for(int k=0; k<ntrunk05; k++) trnksum05 += buffer[k];
for(int k=0; k<ntrunk06; k++) trnksum06 += buffer[k];
for(int k=0; k<ntrunk07; k++) trnksum07 += buffer[k];
for(int k=0; k<ntrunk08; k++) trnksum08 += buffer[k];
for(int k=0; k<ntrunk09; k++) trnksum09 += buffer[k];
for(int k=0; k<ntrunk10; k++) trnksum10 += buffer[k];
ntp[0]=0;
ntp[1]=tecnhits;
ntp[2]=ntrunk06;
ntp[3]=trnksum06;
ntp[4]=dphiproj;
ntp[5]=dphipc3;
ntp[6]=mom;
ntp[7]=toftof;
ntp[8]=tofpl;
ntp[9]=bbct0;
ntp[10]=charge;
ntp[11]=tecalpha;
ntp[12]=dchalpha;
ntp[13]=(float)cent;
ntp[14]=(float)tec->get_TecNTrack();
ntp[15]=(float)cgl->get_CglNTrack();
ntp[16]=trnksum05;
ntp[17]=trnksum07;
ntp[18]=trnksum10;
ntp[19]=ntrunk05;
ntp[20]=ntrunk07;
ntp[21]=ntrunk10;
ntp[22]=mytrklen;
ntp[23]=(float)el;
ntp[24]=(float)npmt0;
ntp[25]=ntrunk07;
ntp[26]=ntrunk08;
ntp[27]=trnksum08;
ntp[28]=trnksum09;
ntp[29]=0.;
ntp[30]=0.;
ntp[31]=1.;
ntp[32]=tec->getTrackNhits(itrk,0);
ntp[33]=tec->getTrackNhits(itrk,1);
ntp[34]=tec->getTrackNhits(itrk,2);
ntp[35]=tec->getTrackNhits(itrk,3);
ntp[36]=tec->getTrackNhits(itrk,4);
ntp[37]=tec->getTrackNhits(itrk,5);
ntp[38]=(float)tecindex;
ntpcheck->Fill(ntp);
// } // good tof
buffer.clear();
} // end loop over tec tracks
HITMAP.clear();
pairs.clear();
evtnum++;
return True;
}
//===========================================================================
PHBoolean mTecPIDModule::analyse(DstContent* dst, TecCalibrationObject* TCO) {
using namespace PhUtilities;
using namespace TecUtilities;
vector<float> buffer;
multimap<int, float, less<int> > HITMAP;
float ntp[39];
// TriggerHelper triggerhelper(dst);
// L2DecisionHelper lvl2Helper(dst);
// TrigLvl1 *trig = static_cast<TrigLvl1*> (dst->getClass("TrigLvl1"));
// TrigRunLvl1* trigrun1 = static_cast<TrigRunLvl1*> (dst->getClass("TrigRunLvl1"));
// TrigRunLvl2* trigrun2 = static_cast<TrigRunLvl2*> (dst->getClass("TrigRunLvl2"));
// Lvl2DecisionOut* lvl2 = static_cast<Lvl2DecisionOut*> (dst->getClass("L2Decision"));
// RunHeader *run = static_cast<RunHeader*> (dst->getClass("RunHeader"));
// EventHeader *evt = static_cast<EventHeader*> (dst->getClass("EventHeader"));
BbcOut *bbc = static_cast<BbcOut*> (dst->getClass("BbcOut"));
// VtxOut *vtx = static_cast<VtxOut*> (dst->getClass("VtxOut"));
ZdcOut *zdc = static_cast<ZdcOut*> (dst->getClass("ZdcOut"));
TecOut *tec = static_cast<TecOut*> (dst->getClass("TecOut"));
DchTrack *dch = static_cast<DchTrack*> (dst->getClass("DchTrack"));
CglTrack* cgl = static_cast<CglTrack*> (dst->getClass("CglTrack"));
// dCglTrackWrapper* cgl = static_cast<dCglTrackWrapper*> (dst->getClass("dCglTrack"));
CrkRing* crk = static_cast<CrkRing*> (dst->getClass("CrkRing"));
// dTofReconstructedWrapper* tof = static_cast<dTofReconstructedWrapper*> (dst->getClass("dTofReconstructed"));
TofOut* tof = static_cast<TofOut*> (dst->getClass("TofOut"));
// PadCluster* pc1 = static_cast<PadCluster*> (dst->getClass("Pc1Cluster"));
// PadCluster* pc2 = static_cast<PadCluster*> (dst->getClass("Pc2Cluster"));
// dPadClusterWrapper* pc3 = static_cast<dPadClusterWrapper*> (dst->getClass("dPc3Cluster"));
PadCluster* pc3 = static_cast<PadCluster*> (dst->getClass("Pc3Cluster"));
// EmcClusterLocalExt* emc =
// static_cast<EmcClusterLocalExt*> (dst->getClass("EmcClusterLocalExt"));
PHTrackOut* proj = static_cast<PHTrackOut*> (dst->getClass("PHTrackOut"));
TecGeometryObject* TGO = 0;
static float beg1=0.;
static float beg2=0.;
if (!init_done) {
cout << "Initializing random number generator..." << endl;
rng = gsl_rng_alloc(gsl_rng_mt19937);
TGO = new TecGeometryObject();
TGO->FetchFromFile();
mTecPIDModule::GetBeg12(TGO, &beg1, &beg2);
delete TGO;
init_done = 1;
}
if(evtnum%100==0) cout << "event # " << evtnum << endl;
// make hit list
for(int j=0; j<tec->getNHits(); j++) {
int trkid=tec->getHitTrackID(j);
if(trkid>-1) {
int jindex = tec->getHitIndex(j);
int myadc = tec->getHitADC(j);
float mycharge = Ampl2Charge(myadc);
float gainA = TCO->getAbsoluteGain(jindex);
mycharge = mycharge*gainA;
HITMAP.insert(pair<int, float>(trkid,mycharge));
}
}
// try to find conversion pairs
vector<int> pairs;
/*
for(int i=0; i<(int)dch->get_DchNTrack()-1; i++) {
for(int j=i+1; j<(int)dch->get_DchNTrack(); j++) {
if(dch->get_arm(i)==0 && dch->get_arm(j)==0) {
if((dch->get_alpha(i)*dch->get_alpha(j)<0)) {
float px1 = dch->get_momentum(i)*sin(dch->get_theta0(i))*cos(dch->get_phi0(i));
float py1 = dch->get_momentum(i)*sin(dch->get_theta0(i))*sin(dch->get_phi0(i));
float pz1 = dch->get_momentum(i)*cos(dch->get_theta0(i));
float px2 = dch->get_momentum(j)*sin(dch->get_theta0(j))*cos(dch->get_phi0(j));
float py2 = dch->get_momentum(j)*sin(dch->get_theta0(j))*sin(dch->get_phi0(j));
float pz2 = dch->get_momentum(j)*cos(dch->get_theta0(j));
float phiv = get_phiv(px1,py1,pz1,px2,py2,pz2);
float invmass = get_invmass(px1,py1,pz1,px2,py2,pz2);
if(fabs(phiv)<0.1 && invmass>0.065 && invmass<0.090) {
pairs.push_back(i); pairs.push_back(j);
}
}
}
}
}
*/
float bbct0 = bbc->get_TimeZero();
float bbc1 = bbc->get_ChargeSum(Bbc::South);
float bbc2 = bbc->get_ChargeSum(Bbc::North);
float zdc1 = zdc->get_Energy(Zdc::South);
float zdc2 = zdc->get_Energy(Zdc::North);
int cent = getCentralityByClock(bbc1, bbc2, zdc1, zdc2);
for(int itrk=0; itrk<(int)tec->get_TecNTrack(); itrk++) {
float dchalpha = -999.;
float tecalpha = tec->getTrackAlpha(itrk);
float tecphi = tec->getTrackPhi(itrk);
float tecnhits = (float)tec->getTrackNhits(itrk);
float pc3phi = 9999.;
float projphi = 9999.;
float dphiproj = 9999.;
float dphipc3 = 9999.;
float mom = 0.;
float toftof = 9999.;
float tofpl = 999.;
float charge = 1.;
float theta0 = 0.;
// float trklen = tec->getTrackLength(itrk);
float mytrklen = -999.;
int el = 0;
int npmt0 = 0;
float pc3z=999.;
float projz=999.;
float zcorr=999.;
int tecindex = tec->getTrackIndex(itrk);
// calculate track length
float dx = tec->getTrackXout(itrk)-tec->getTrackXin(itrk);
float dy = tec->getTrackYout(itrk)-tec->getTrackYin(itrk);
float dz = 0.0;
// Check if this Tec track was associated with a global track
// and calculate phi difference for Dch track projection and Pc3 cluster
for(int icgl=0; icgl<(int)cgl->get_CglNTrack(); icgl++) {
if(cgl->get_tectrackid(icgl)==itrk) {
float xp = proj->get_projectionTec(icgl,0);
float yp = proj->get_projectionTec(icgl,1);
projz = proj->get_projectionTec(icgl,2);
if(xp!=0) projphi = 3.1415927+atan(yp/xp);
if(projphi<0) projphi=2*3.1415927+projphi;
dphiproj = tecphi - projphi;
int pc3ptr = cgl->get_pc3clusid(icgl);
if(pc3ptr>-1 && pc3ptr<(int)pc3->get_PadNCluster()) {
float x = pc3->get_xyz(0,pc3ptr);
float y = pc3->get_xyz(1,pc3ptr);
pc3z = pc3->get_xyz(2,pc3ptr);
if(x!=0) pc3phi = 3.1415927+atan(y/x);
if(pc3phi<0) pc3phi=2*3.1415927+pc3phi;
dphipc3 = tecphi-pc3phi;
}
if(dch->get_alpha(icgl)>0) charge = -1.;
mom = dch->get_momentum(icgl);
dchalpha = dch->get_alpha(icgl);
theta0 = dch->get_theta0(icgl);
dz = fabs(beg1-beg2)*tan(3.1415927/2.-theta0);
int tofptr = cgl->get_tofrecid(icgl);
if(tofptr>-1) {
toftof = tof->get_tof(tofptr);
tofpl = proj->get_tofPathLength(icgl);
}
int crkptr = cgl->get_richringid(icgl);
if(crkptr>-1) {
npmt0 = crk->get_npmt0(crkptr);
}
// // check if this is a conversion electron
// for(int i=0; i<(int)pairs.size(); i++) {
// if(icgl==pairs[i] && npmt0>0) el = 1;
// }
break;
}
}
mytrklen = sqrt(dx*dx+dy*dy+dz*dz);
// if(toftof<999.) {
//cout << "trklen: " << trklen << " " << mytrklen << " "
// << dx << " " << dy << " " << dz << " "
// << 180./3.1415927*(3.1415927/2.-theta0) << " "
// << projz << " " << beg1-beg2 << " " << mom << endl;
// if(pc3z!=999.) {
// zcorr = 9.03343e-01+2.02908e-03*pc3z-1.05649e-05*pc3z*pc3z;
// }
// else {
if(projz!=999.) {
zcorr = 9.03343e-01+2.02908e-03*fabs(projz)-1.05649e-05*projz*projz;
}
else {zcorr=1.;}
// }
// Calculate truncated dE/dX
// Fill buffer from hitmap
for(multimap<int, float, less<int> >::iterator it = HITMAP.begin(); it != HITMAP.end(); ++it)
{
if((*it).first==itrk) {
buffer.push_back((*it).second);
}
}
// Re-order buffer
// Re-order buffer
sort(buffer.begin(), buffer.end());
int ntrunk05 = (int)(buffer.size()*0.5);
int ntrunk06 = (int)(buffer.size()*0.6);
int ntrunk07 = (int)(buffer.size()*0.7);
int ntrunk08 = (int)(buffer.size()*0.8);
int ntrunk09 = (int)(buffer.size()*0.9);
int ntrunk10 = (int)buffer.size();
float trnksum05 = 0.;
float trnksum06 = 0.;
float trnksum07 = 0.;
float trnksum08 = 0.;
float trnksum09 = 0.;
float trnksum10 = 0.;
for(int k=0; k<ntrunk05; k++) trnksum05 += buffer[k];
for(int k=0; k<ntrunk06; k++) trnksum06 += buffer[k];
for(int k=0; k<ntrunk07; k++) trnksum07 += buffer[k];
for(int k=0; k<ntrunk08; k++) trnksum08 += buffer[k];
for(int k=0; k<ntrunk09; k++) trnksum09 += buffer[k];
for(int k=0; k<ntrunk10; k++) trnksum10 += buffer[k];
ntp[0]=(float)evtnum;
ntp[1]=(float)tecnhits;
ntp[2]=ntrunk06;
ntp[3]=trnksum06;
ntp[4]=dphiproj;
ntp[5]=dphipc3;
ntp[6]=mom;
ntp[7]=toftof;
ntp[8]=tofpl;
ntp[9]=bbct0;
ntp[10]=charge;
ntp[11]=tecalpha;
ntp[12]=dchalpha;
ntp[13]=(float)cent;
ntp[14]=(float)tec->get_TecNTrack();
ntp[15]=(float)cgl->get_CglNTrack();
ntp[16]=trnksum05;
ntp[17]=trnksum07;
ntp[18]=trnksum10;
ntp[19]=ntrunk05;
ntp[20]=ntrunk07;
ntp[21]=ntrunk10;
ntp[22]=mytrklen;
ntp[23]=(float)el;
ntp[24]=(float)npmt0;
ntp[25]=ntrunk07;
ntp[26]=ntrunk08;
ntp[27]=trnksum08;
ntp[28]=trnksum09;
ntp[29]=pc3z;
ntp[30]=projz;
ntp[31]=zcorr;
ntp[32]=tec->getTrackNhits(itrk,0);
ntp[33]=tec->getTrackNhits(itrk,1);
ntp[34]=tec->getTrackNhits(itrk,2);
ntp[35]=tec->getTrackNhits(itrk,3);
ntp[36]=tec->getTrackNhits(itrk,4);
ntp[37]=tec->getTrackNhits(itrk,5);
ntp[38]=(float)tecindex;
ntpcheck->Fill(ntp);
cout << "good tof: " << toftof << " " << tofpl << " " << bbct0
<< " " << mom << " " << trnksum05 << endl;
// } // good tof
buffer.clear();
} // end loop over tec tracks
HITMAP.clear();
pairs.clear();
evtnum++;
return True;
}
//===========================================================================
PHBoolean mTecPIDModule::analyse(TChain* chain1, TChain* chain2, TecCalibrationObject* TCO) {
using namespace PhUtilities;
using namespace TecUtilities;
//chain1->Print();
TFile* OutputNtupleFile = new TFile("test.root","RECREATE");
TNtuple* ntpcheck = new TNtuple("ntpcheck"," ","evt:nhits:ntrunk:trnksum:dphiproj:dphipc3:mom:toftof:tofpl:bbct0:charge:tecalpha:dchalpha:cent:tecmult:dchmult:trnksum05:trnksum07:trnksum10:ntrunk05:ntrunk07:ntrunk10");
float ntp[22];
vector<float> buffer;
multimap<int, float, less<int> > HITMAP;
rng = gsl_rng_alloc(gsl_rng_mt19937);
// EventHeader* evt = 0;
// RunHeader* run = 0;
// PadCluster* pc3 = 0;
dPadClusterWrapper* pc3 = 0;
// CglTrack* cgl = 0;
dCglTrackWrapper* cgl = 0;
dTofReconstructedWrapper* tof = 0;
DchTrack* dch = 0;
TecOut* tec = 0;
PHTrackOut* proj = 0;
BbcOut* bbc = 0;
ZdcOut* zdc = 0;
cout << "setting branch addresses..." << endl;
// chain2->SetBranchAddress("RUN/RunHeader", &run);
// chain1->SetBranchAddress("DST/EventHeader", &evt);
chain1->SetBranchAddress("DST/dPc3Cluster", &pc3);
chain1->SetBranchAddress("DST/DchTrack", &dch);
chain1->SetBranchAddress("DST/dCglTrack", &cgl);
chain1->SetBranchAddress("DST/dTofReconstructed", &tof);
chain1->SetBranchAddress("DST/TecOut", &tec);
chain1->SetBranchAddress("DST/PHTrackOut", &proj);
chain1->SetBranchAddress("DST/BbcOut", &bbc);
chain1->SetBranchAddress("DST/ZdcOut", &zdc);
// run info
// chain2->GetEvent(1);
cout << "Total events: " << chain1->GetEntries() << endl;
// loop over events
for(int i=1; i<chain1->GetEntries(); i++) {
if(i%1==0) cout << "event # " << i << endl;
cout << "***getting event..." << endl;
chain1->GetEvent(i);
cout << "***done." << endl;
// make hit list
for(int j=0; j<tec->getNHits(); j++) {
int trkid=tec->getHitTrackID(j);
if(trkid>-1) {
int jindex = tec->getHitIndex(j);
int myadc = tec->getHitADC(j);
float mycharge = Ampl2Charge(myadc);
float gainA = TCO->getAbsoluteGain(jindex);
mycharge = mycharge*gainA;
HITMAP.insert(pair<int, float>(trkid,mycharge));
}
}
cout << "***centrality..." << endl;
float bbct0 = bbc->get_TimeZero();
float bbc1 = bbc->get_ChargeSum(Bbc::South);
float bbc2 = bbc->get_ChargeSum(Bbc::North);
float zdc1 = zdc->get_Energy(Zdc::South);
float zdc2 = zdc->get_Energy(Zdc::North);
int cent = getCentralityByClock(bbc1, bbc2, zdc1, zdc2);
for(int itrk=0; itrk<(int)tec->get_TecNTrack(); itrk++) {
float dchalpha = -999.;
float tecalpha = tec->getTrackAlpha(itrk);
float tecphi = tec->getTrackPhi(itrk);
float tecnhits = (float)tec->getTrackNhits(itrk);
// int tecsect = tec->getTrackSector(itrk);
float pc3phi = 9999.;
float projphi = 9999.;
float dphiproj = 9999.;
float dphipc3 = 9999.;
float mom = 0.;
float toftof = 9999.;
float tofpl = 999.;
float charge = 1.;
// Check if this Tec track was associated with a global track
// and calculate phi difference for Dch track projection and Pc3 cluster
for(int icgl=0; icgl<(int)cgl->RowCount(); icgl++) {
if(cgl->get_tectrackid(icgl)==itrk) {
float xp = proj->get_projectionTec(icgl,0);
float yp = proj->get_projectionTec(icgl,1);
if(xp!=0) projphi = 3.1415927+atan(yp/xp);
if(projphi<0) projphi=2*3.1415927+projphi;
dphiproj = tecphi - projphi;
int pc3ptr = cgl->get_pc3clusid(icgl);
if(pc3ptr>-1 && pc3ptr<(int)pc3->RowCount()) {
float x = pc3->get_xyz(0,pc3ptr);
float y = pc3->get_xyz(1,pc3ptr);
if(x!=0) pc3phi = 3.1415927+atan(y/x);
if(pc3phi<0) pc3phi=2*3.1415927+pc3phi;
dphipc3 = tecphi-pc3phi;
}
if(dch->get_alpha(icgl)>0) charge = -1.;
mom = dch->get_momentum(icgl);
dchalpha = dch->get_alpha(icgl);
int tofptr = cgl->get_tofrecid(icgl);
if(tofptr>-1) {
toftof = tof->get_tof(tofptr);
}
tofpl = proj->get_tofPathLength(icgl);
break;
}
}
if(toftof<999.) {
// Calculate truncated dE/dX
// Fill buffer from hitmap
for(multimap<int, float, less<int> >::iterator it = HITMAP.begin(); it != HITMAP.end(); ++it)
{
if((*it).first==itrk) {
buffer.push_back((*it).second);
}
}
// Re-order buffer
sort(buffer.begin(), buffer.end());
int ntrunk05 = (int)(buffer.size()*0.5);
int ntrunk06 = (int)(buffer.size()*0.6);
int ntrunk07 = (int)(buffer.size()*0.7);
int ntrunk10 = buffer.size();
float trnksum05 = 0.;
float trnksum06 = 0.;
float trnksum07 = 0.;
float trnksum10 = 0.;
for(int k=0; k<ntrunk05; k++) trnksum05 += buffer[k];
for(int k=0; k<ntrunk06; k++) trnksum06 += buffer[k];
for(int k=0; k<ntrunk07; k++) trnksum07 += buffer[k];
for(int k=0; k<ntrunk10; k++) trnksum10 += buffer[k];
ntp[0]=(float)i;
ntp[1]=(float)tecnhits;
ntp[2]=ntrunk06;
ntp[3]=trnksum06;
ntp[4]=dphiproj;
ntp[5]=dphipc3;
ntp[6]=mom;
ntp[7]=toftof;
ntp[8]=tofpl;
ntp[9]=bbct0;
ntp[10]=charge;
ntp[11]=tecalpha;
ntp[12]=dchalpha;
ntp[13]=(float)cent;
ntp[14]=(float)tec->get_TecNTrack();
ntp[15]=(float)cgl->RowCount();
ntp[16]=trnksum05;
ntp[17]=trnksum07;
ntp[18]=trnksum10;
ntp[19]=ntrunk05;
ntp[20]=ntrunk07;
ntp[21]=ntrunk10;
ntpcheck->Fill(ntp);
} // good tof
buffer.clear();
} // end loop over tec tracks
HITMAP.clear();
} // end loop over events
OutputNtupleFile->Write();
OutputNtupleFile->Close();
return True;
}
//===========================================================================
PHBoolean mTecPIDModule::event(PHCompositeNode *root,
TecCalibrationObject* TCO) {
if(Verbose>0) cout << "mTecPIDModule: Started..." << endl;
PHNodeIterator iii(root);
//------ Get pointers to the Tables ------------------------------------
dDchTracksWrapper* dDchTracks = 0;
dDchTracksNode_t* DCHTN = static_cast<dDchTracksNode_t*>(iii.findFirst("PHIODataNode", "dDchTracks"));
if (!DCHTN) {
if(Verbose>0)
cerr << "mTecPIDModule -> WARNING: dDchTracks table not found." << endl;
}
else {
dDchTracks = DCHTN->getData();
}
dCglTrackWrapper* dCglTrack = 0;
dCglTrackNode_t* CGLTN = static_cast<dCglTrackNode_t*>(iii.findFirst("PHIODataNode", "dCglTrack"));
if (!CGLTN) {
if(Verbose>0)
cerr << "mTecPIDModule -> WARNING: dCglTrack table not found." << endl;
}
else {
dCglTrack = CGLTN->getData();
}
//-------------------------------------------------------------------------------
// TecOutV1
PHTypedNodeIterator<TecOutV1> teciter(root);
TecOutNode_t *TecOutNode = teciter.find("TecOutV1");
TecOutV1* tecout;
if(TecOutNode) {
tecout = (TecOutV1*)TecOutNode->getData();
}
else {
if(Verbose>0) cerr << "mTecPIDModule ERROR: Can not find TecOutV1." << endl;
return False;
}
// TecOutV2
PHTypedNodeIterator<TecOut> tecitershort(root);
TecOutNodeShort_t *TecOutNodeShort = tecitershort.find("TecOut");
TecOut* tecoutshort;
if(TecOutNodeShort) {
tecoutshort = (TecOutV2*)TecOutNodeShort->getData();
}
else {
if(Verbose>0) cerr << "mTecPIDModule ERROR: Can not find TecOut." << endl;
return False;
}
//--------------------------------------------------------------------
static int EventNumber = 0;
vector<float> buffer;
multimap<int, float, less<int> > HITMAP;
if(Verbose > 0) {
cout << "mTecPID -> Started event # " << EventNumber << endl;
cout << "mTecPID -> Number of tracks: " << tecout->getNTracks() << " "
<< tecoutshort->getNTracks() << endl;
cout << "mTecPID -> Number of hits: " << tecout->getNHits() << endl;
if(dCglTrack) printf("mTecPID -> # of entries in dCglTrack: %d\n",dCglTrack->RowCount());
if(dDchTracks) printf("mTecPID -> # of entries in dDchTracks: %d\n",dDchTracks->RowCount());
}
// Fill HITMAP only if number of hits in a plane is OK
int ntmp=0;
for(int j=0; j<tecout->getNHits(); j++) {
//int jbin = tecout->getHitTimeBin(j);
//int jindex = tecout->getHitIndex(j);
float myampl = tecout->getHitCharge(j);
//float gainT = TCO->getTimingGain(jindex,jbin);
//myampl = myampl*gainT;
tecout->setHitCharge(j,myampl);
int trkid=tecout->getHitTrackID(j);
if(trkid>-1) {
if(tecout->getTrackNhits(trkid,tecout->getHitPlane(j))>hitcut &&
tecout->getTrackNhits(trkid,tecout->getHitPlane(j))<80) {
ntmp++;
HITMAP.insert(pair<int, float>(trkid,myampl));
}
}
} // j - end loop over hits
if(Verbose > 0) cout << "mTecPID -> Number of associated hits = " << ntmp << endl;
//======================= Loop over TEC tracks =========================
for(int itrk=0; itrk<tecout->getNTracks(); itrk++) {
// Find track length
float TrackLength = 0.;
if(dCglTrack && dDchTracks) {
int status = GetTrackLength(itrk, tecout, dCglTrack, dDchTracks, TrackLength);
}
// Find number of fired planes for this track
int nfpl=0;
for(int ipl=0; ipl<TECMAXPLANE; ipl++) {
if(tecout->getTrackNhits(itrk,ipl)>hitcut &&
tecout->getTrackNhits(itrk,ipl)<80) {nfpl++;}
}
// Correct TrackLength for number of fired planes
TrackLength = TrackLength/float(TECMAXPLANE)*float(nfpl);
// Find dE/dX for this track
// Fill buffer first
for(multimap<int, float, less<int> >::iterator it = HITMAP.begin();
it != HITMAP.end(); ++it) {
if((*it).first==itrk) {
buffer.push_back((*it).second);
}
}
if(buffer.size()>10) {
// Re-order buffer
sort(buffer.begin(), buffer.end());
// Calculate truncated mean dE/dX
int ntrunk = ((int)(buffer.size()*Truncation));
float trnksum = 0.; float trnksumtot=0.;
for(int i=0; i<ntrunk; i++) trnksum += buffer[i];
for(int i=0; i<(int)buffer.size(); i++) trnksumtot += buffer[i];
float trnksum_perbin=0.;
float trnksumtot_perbin=0.;
float trnksum_perdx=0.;
float trnksumtot_perdx=0.;
if(ntrunk > 0) {
trnksum_perbin = trnksum / ntrunk;
}
trnksumtot_perbin = trnksumtot / buffer.size();
if(TrackLength!=0.) {
trnksum_perdx=trnksum/TrackLength;
trnksumtot_perdx=trnksumtot/TrackLength;
}
// Fill TecOut
tecout->setTrackdEdX(itrk, trnksum);
tecoutshort->setTrackdEdX(itrk, trnksum);
tecout->setTrackNdEdXbins(itrk, ntrunk);
tecoutshort->setTrackNdEdXbins(itrk, ntrunk);
tecout->setTrackLength(itrk, TrackLength);
tecoutshort->setTrackLength(itrk, TrackLength);
if(Verbose>1) {
cout << " track: " << itrk << endl;
cout << " nbins: " << ntrunk << endl;
cout << " nbinstot: " << buffer.size() << endl;
cout << " FADC_sum: " << tecoutshort->getTrackdEdX(itrk) << endl;
cout << " FADC_sumtot: " << trnksumtot << endl;
cout << " Truncation: " << Truncation << endl;
cout << " TrackLength: " << tecoutshort->getTrackLength(itrk) << endl;
}
buffer.clear(); // reset buffer
} // buffer sixe > 10
} //========================= itrk - end loop over TEC tracks
EventNumber++;
HITMAP.clear();
if(Verbose>0) { cout << "mTecPIDModule: Finished." << endl; }
return True;
}
//===========================================================================
PHBoolean mTecPIDModule::event(PHCompositeNode *root) {
PHNodeIterator iii(root);
if(Verbose>0) cout << "mTecPIDModule: Started..." << endl;
if(Verbose>1) cout << "mTecPIDModule -> Instantiating PhHistogramFactory..." << endl;
PhHistogramFactory *hf = PhHistogramFactory::instance();
static TNtuple *ntp0;
float ntpart[8];
int foundDch=1; int foundCglT=1; int foundCglP=1;
static int first = 0;
if(first==0) {
// Book the Ntuple
first=1;
hf->cd();
ntp0 = new TNtuple("ntp0","tecpid",
"evt:ptot:nbins:trklen:dedx:gpart:tecmom:teccharge");
}
else
{
hf->cd();
}
//------ Get pointers to the Tables ------------------------------------
dDchTracksNode_t* DCHTN = static_cast<dDchTracksNode_t*>(iii.findFirst("PHIODataNode", "dDchTracks"));
if (!DCHTN) {
if(Verbose>0)
cerr << "mTecPIDModule -> WARNING: dDchTracks table not found." << endl;
foundDch=0;
}
dDchTracksWrapper* dDchTracks = 0;
if(foundDch) dDchTracks = DCHTN->getData();
dCglTrackNode_t* CGLTN = static_cast<dCglTrackNode_t*>(iii.findFirst("PHIODataNode", "dCglTrack"));
if (!CGLTN) {
if(Verbose>0)
cerr << "mTecPIDModule -> WARNING: dCglTrack table not found." << endl;
foundCglT=0;
}
dCglTrackWrapper* dCglTrack = 0;
if(foundCglT) dCglTrack = CGLTN->getData();
dCglParticleNode_t* CGLPN = static_cast<dCglParticleNode_t*>(iii.findFirst("PHIODataNode", "dCglParticle"));
if (!CGLPN) {
if(Verbose>0)
cerr << "mTecPIDModule -> WARNING: dCglParticle table not found." << endl;
foundCglP=0;
}
dCglParticleWrapper* dCglParticle = 0;
if(foundCglP) dCglParticle = CGLPN->getData();
//-------------------------------------------------------------------------------
// TecOutV1
PHTypedNodeIterator<TecOutV1> teciter(root);
TecOutNode_t *TecOutNode = teciter.find("TecOutV1");
TecOutV1* tecout;
if(TecOutNode) {
tecout = (TecOutV1*)TecOutNode->getData();
}
else {
if(Verbose>0) cerr << "mTecPIDModule ERROR: Can not find TecOutV1." << endl;
return False;
}
// TecOutV2
PHTypedNodeIterator<TecOut> tecitershort(root);
TecOutNodeShort_t *TecOutNodeShort = tecitershort.find("TecOut");
TecOut* tecoutshort;
if(TecOutNodeShort) {
tecoutshort = (TecOutV2*)TecOutNodeShort->getData();
}
else {
if(Verbose>0) cerr << "mTecPIDModule ERROR: Can not find TecOut." << endl;
return False;
}
// Find TecCalibrationObject and TecGeometryObject in the node tree
PHBoolean status2 = False;
PHBoolean status3 = False;
TecCalibrationObject *TCO;
TecGeometryObject *TGO;
ObjectNode_t *TCONode;
ObjectNode_t *TGONode;
PHNodeIterator nodeIter (root), *parNodeIter;
PHCompositeNode *parNode;
parNode =
static_cast < PHCompositeNode * >(nodeIter.findFirst ("PHCompositeNode", "PAR"));
if (!parNode) {
PHMessage ("mTecPIDModule: ", PHError,
"PAR node does not exist.");
return False;
}
parNodeIter = new PHNodeIterator (parNode);
TCONode =
static_cast < ObjectNode_t * >(parNodeIter->findFirst ("PHIODataNode", "TecCalibration"));
if (!TCONode) {
PHMessage ("mTecPIDModule: ", PHError,
"TecCalibrationObject not found.\n");
PHMessage ("mTecPIDModule: ", PHError,
"Will set all Calibration Constatnts to 1.\n");
TCO = new TecCalibrationObject ();
status2 = TCO->FetchFromFile();
}
else {
TCO = static_cast < TecCalibrationObject * >(TCONode->getData ());
}
TGONode =
static_cast < ObjectNode_t * >(parNodeIter->findFirst ("PHIODataNode", "TecGeometry"));
if (!TGONode) {
PHMessage ("mTecPIDModule: ", PHError,
"TecGeometryObject not found.\n");
PHMessage ("mTecPIDModule: ", PHError,
"Will set Tec Geometry to default Geant Geometry.\n");
TGO = new TecGeometryObject ();
status3 = TGO->FetchFromFile();
}
else {
TGO = static_cast < TecGeometryObject * >(TGONode->getData ());
}
delete parNodeIter;
//--------------------------------------------------------------------
int Debug,itrk,dchtrkid,ntrunk,status = 0,nbtgm,mypart = 0;
int nok,Particle;
float trnksum,trnksumtot;
float Momentum,TrackLength,betagamma = 0.0;
float rtrnksum,rtrnksumtot,x,m,Gain,calibration,truncation;
float gaincorr,gaincorrel,trnkcorr,trnkcorrel,chanperbin;
float e_separation,prob_e,prob_pi,prob_p,prob_k,btgm[2],min;
float mean1,mean2,mean3,mean4,sigma1,sigma2,sigma3,sigma4;
float tmprtrnksum = 0.0;
FILE *out_file = 0;
static int EventNumber = 0;
vector<float> buffer; // stl vector of floats
int ibuf,ibuf0;
multimap<int, float, less<int> > HITMAP;
Debug = Verbose;
// Open ASCII file for writing results
if (Write2File == 1)
{
out_file = fopen("pid.dat","w");
}
if(Debug > 0) {
cout << "mTecPID -> Started event # " << EventNumber << endl;
cout << "mTecPID -> Number of tracks: " << tecout->getNTracks() << " "
<< tecoutshort->getNTracks() << endl;
cout << "mTecPID -> Number of hits: " << tecout->getNHits() << endl;
if(foundCglT) printf("mTecPID -> # of entries in dCglTrack: %d\n",dCglTrack->RowCount());
if(foundCglP) printf("mTecPID -> # of entries in dCglParticle: %d\n",dCglParticle->RowCount());
}
// Make a list of geant tracks if perfect tracking was requested
int gtrknum0 = 0;
if (PerfectPID != 0)
{
gtrknum0 = MakeGTrackList(root);
}
int ntmp=0;
// Apply Timing Calibration and fill HITMAP
for(int j=0; j<tecout->getNHits(); j++) {
int jbin = tecout->getHitTimeBin(j);
int jindex = tecout->getHitIndex(j);
float myampl = tecout->getHitCharge(j);
float gainT = TCO->getTimingGain(jindex,jbin);
myampl = myampl*gainT;
tecout->setHitCharge(j,myampl);
int trkid=tecout->getHitTrackID(j);
if(trkid>-1) {
ntmp++;
HITMAP.insert(pair<int, float>(trkid,myampl));
}
} // j - end loop over hits
if(Debug > 0) cout << "mTecPID -> Number of associated hits = " << ntmp << endl;
// Loop over TEC tracks
nok=0;
for(itrk=0; itrk<tecout->getNTracks(); itrk++) {
ibuf=0;
ibuf0=0;
// Find transverse momentum and charge from alpha angle in Tec
float tecPt, tecCharge, tecMomentum;
GuessMomentum(itrk, root, &tecPt, &tecCharge);
// Try to find total Tec momentum using Z from Pc1 and Pc3
int pc1id=-1;
int pc3id=-1;
float DZ=-999.;
int statdz;
statdz = CalculateDZ(root, itrk, Debug, &pc1id, &pc3id, &DZ);
if(statdz==0 && DZ!=-999.) {
tecMomentum = tecPt*sqrt(1.+DZ*DZ);
}
else {
tecMomentum = tecPt;
}
// Is there a global track associated with this Tec track?
// If yes, find Momentum from global tracking.
Momentum = -1.;
Particle = -1;
if(PerfectPID==0) {
if(foundCglT && foundCglP) {
status = GetMomentum(itrk, Debug, dCglTrack, dCglParticle,
&Momentum, &dchtrkid);
}
}
else {
status = GetPerfectMomentum(itrk, root, gtrknum0,
&Momentum, &Particle, Debug);
}
// If using of Tec Pt is requested, or no global track is associated, use tecMomentum
if(UseTec>0 || Momentum<0.) { Momentum = tecMomentum; }
if(Debug>1) printf("mTecPID -> after GetMomentum: %d %f\n",status,Momentum);
if(Momentum>0.) {
if(UseTec>0 && PerfectPID==0) method = 1; // Tec-only information can not give track length
// Use dE/dX per time bin;
// Find track length
TrackLength = -999.;
int pc1id=-1;
int pc3id=-1;
if(PerfectPID==0) {
status = GetTrackLengthFromPAD(itrk, root, Debug,
&pc1id, &pc3id, &TrackLength);
if(TrackLength==-999.) {
if(foundCglT && foundCglP && foundDch) {
status = GetTrackLength(itrk, Debug, tecout, TGO,
dCglTrack, dDchTracks,
&dchtrkid, &TrackLength);
}
}
}
else {
status = GetPerfectTrackLength(itrk, root, gtrknum0,
&TrackLength, Debug);
}
if(Debug>1) cout << "mTecPID -> Track Length found." << endl;
// Find dE/dX for this track
// Fill buffer
for(multimap<int, float, less<int> >::iterator it = HITMAP.begin();
it != HITMAP.end(); ++it) {
if((*it).first==itrk) {
buffer.push_back((*it).second);
}
}
ibuf=buffer.size();
ibuf0=buffer.size();
// Re-order buffer
sort(buffer.begin(), buffer.end());
// Print buffer contents
if(Debug > 1) {
cout << "buffer size = "<<buffer.size()<<" "<<ibuf<<" "<<ibuf0<<endl;
}
if(Debug > 20) {
for (int i = 0; i < (int)buffer.size(); i++) {
cout << buffer[i] << " ";
}
cout << endl;
}
// Calculate truncated mean dE/dX
truncation = Truncation;
ntrunk = ((int)(ibuf*truncation));
trnksum = 0.; trnksumtot=0.;
for(int i=0; i<ntrunk; i++) trnksum += buffer[i];
for(int i=0; i<ibuf; i++) trnksumtot += buffer[i];
rtrnksum=0.;
rtrnksumtot=0.;
if (ntrunk > 0)
{
tmprtrnksum = trnksum / ntrunk;
}
// use dE/dX per one time bin
if(method==1) {
if(ntrunk>0) {
rtrnksum=trnksum/ntrunk;
rtrnksumtot=trnksumtot/ibuf;
}
else {
rtrnksum=999.;
rtrnksumtot=999.;
}
}
// use dE/dX per unit of track length
// ibuf0 is the number of bins without rejecting overlaps, ibuf is the same
// number when overlaps are rejected. If RejectOverlaps==0, they are equal.
// 0.65 is empirical factor, needs to be checked.
if(method==0) {
if((TrackLength!=0) && (TrackLength!=-999.) && (ibuf!=0)) {
rtrnksum=((float)(ibuf0))/((float)(ibuf))*trnksum/TrackLength*0.65;
rtrnksumtot=((float)(ibuf0))/((float)(ibuf))*trnksumtot/TrackLength*0.65;
}
else {
rtrnksum=999.;
rtrnksumtot=999.;
}
}
if(Debug > 1) printf("mTecPID -> ****** TRUNK: %d %d %f %f %f %f\n",itrk,ntrunk,trnksum,trnksumtot,rtrnksum,rtrnksumtot);
if(Debug > 1 && ntrunk > 0) printf("mTecPID -> ****** %f %f \n",trnksum/ntrunk,rtrnksum);
// Calculate probabilities for various particles
calibration = Calibration[0];
Gain = Calibration[1]/1000.;
gaincorr=mip_gaindep[0]+
mip_gaindep[1]*Gain+
mip_gaindep[2]*Gain*Gain;
trnkcorr=mip_trnkdep[0]+
mip_trnkdep[1]*truncation+
mip_trnkdep[2]*truncation*truncation;
gaincorrel=(el_gaindep[0]+
el_gaindep[1]*Gain+
el_gaindep[2]*Gain*Gain)/
(el_gaindep[0]+
el_gaindep[1]*10.+
el_gaindep[2]*100.);
trnkcorrel=el_trnkdep[0]+
el_trnkdep[1]*truncation+
el_trnkdep[2]*truncation*truncation;
chanperbin=mip*gaincorr*trnkcorr*calibration;
e_separation=gaincorrel*trnkcorrel;
if(Debug > 1) printf("mTecPID -> GasGain,etc. = %f %f %f %f %f %f %f\n",
Gain*1000,chanperbin,e_separation,gaincorr,trnkcorr,gaincorrel,trnkcorrel);
// mean* -> expected dE/dX for various particles, rtrnksum -> actual dE/dX
x=Momentum;
mean1=e_separation*chanperbin;
sigma1=mean1*Sigma[0];
prob_e = (rtrnksum-mean1)/sigma1;
m=0.1396;
mean2 = chanperbin*0.0705*(1+(m/x)*(m/x))*
(11.45+log((x/m)*(x/m))-(x/m)*(x/m)/(1+(x/m)*(x/m)));
sigma2 = mean2*Sigma[1];
prob_pi = (rtrnksum-mean2)/sigma2;
m=0.9383;
mean3 = chanperbin*0.0705*(1+(m/x)*(m/x))*
(11.45+log((x/m)*(x/m))-(x/m)*(x/m)/(1+(x/m)*(x/m)));
sigma3 = mean3*Sigma[2];
prob_p = (rtrnksum-mean3)/sigma3;
m=0.4937;
mean4 = chanperbin*0.0705*(1+(m/x)*(m/x))*
(11.45+log((x/m)*(x/m))-(x/m)*(x/m)/(1+(x/m)*(x/m)));
sigma4 = mean4*Sigma[3];
prob_k = (rtrnksum-mean4)/sigma4;
if(Debug > 1) {
printf("mTecPID -> mean_(e,pi,p,k) = %f %f %f %f\n",
mean1,mean2,mean3,mean4);
printf("mTecPID -> prob_(e,pi,p,k) = %f %f %f %f \n",
prob_e,prob_pi,prob_p,prob_k);
}
// Find most probable particle -> particle with smallest dE/dX deviation
min = 9999.;
if(((fabs)(prob_k))<min) {
min = (fabs)(prob_k); mypart = 12;
betagamma = Momentum/0.4937;
}
if(((fabs)(prob_p))<min) {
min = (fabs)(prob_p); mypart = 14;
betagamma = Momentum/0.9383;
}
if(((fabs)(prob_e))<min) {
min = (fabs)(prob_e); mypart = 3;
betagamma = Momentum/0.000511;
}
if(((fabs)(prob_pi))<min) {
min = (fabs)(prob_pi); mypart = 9;
betagamma = Momentum/0.1396;
}
if(Debug > 1) printf("mTecPID -> PARTICLE = %f %d %f %f \n",
rtrnksum,mypart,Momentum,betagamma);
nbtgm=0;
btgm[0]=0.;
btgm[1]=0.;
status = GetBetaGamma(chanperbin, rtrnksum, &btgm[0], &nbtgm);
if(nbtgm==0) {
nbtgm=1;
btgm[0]=3.6;
btgm[1]=0.0;
}
// Write results to a file for calibration
if (Write2File==1 && out_file != 0)
{
status = fprintf(out_file," %f %f %f \n",rtrnksum,tmprtrnksum,Momentum);
}
// Fill TecOut
tecout->setTrackdEdX(itrk, trnksum);
tecoutshort->setTrackdEdX(itrk, trnksum);
tecout->setTrackNdEdXbins(itrk, ntrunk);
tecoutshort->setTrackNdEdXbins(itrk, ntrunk);
tecout->setTrackLength(itrk, TrackLength);
tecoutshort->setTrackLength(itrk, TrackLength);
if(Verbose>1) {
cout << "track: " << itrk << endl;
cout << "nbins: " << ntrunk << endl;
cout << "nbinstot: " << ibuf << endl;
cout << "FADC_sum: " << tecoutshort->getTrackdEdX(itrk) << endl;
cout << "FADC_sumtot: " << trnksumtot << endl;
cout << "Momentum: " << Momentum << endl;
cout << "tecMomentum: " << tecMomentum << endl;
cout << "tecPt: " << tecPt << endl;
cout << "DZ: " << DZ << endl;
cout << "tecCharge: " << tecCharge << endl;
cout << "Truncation: " << Truncation << endl;
cout << "betagamma0: " << btgm[0] << endl;
cout << "betagamma1: " << btgm[1] << endl;
cout << "TrackLength: " << tecoutshort->getTrackLength(itrk) << endl;
}
nok++;
// Fill the Ntuple
if(Write2Ntuple!=0) {
ntpart[0]=EventNumber;
ntpart[1]=Momentum;
ntpart[2]=ntrunk;
ntpart[3]=TrackLength;
ntpart[4]=trnksum;
ntpart[5]=Particle;
ntpart[6]=tecMomentum;
ntpart[7]=tecCharge;
ntp0->Fill(ntpart);
}
buffer.clear(); // reset buffer
} // Momentum > 0
} // itrk - end loop over TEC tracks
if (Write2File == 1 && out_file != 0)
{
status = fclose(out_file);
}
// Save the Ntuple
if(Write2Ntuple!=0) {
hf->cd();
hf->Save("tecpid.root");
}
EventNumber++;
if(Verbose>0) {
cout << "mTecPIDModule: Finished." << endl;
}
// free memory
if (status2) delete TCO;
if (status3) delete TGO;
HITMAP.clear();
return True;
}