//////////////////////////////////////////////////////////////////////////
// //
// ATLFast TrackMaker class. //
// //
// No Pion smearing yet. Muon parameterisations used //
// //
//////////////////////////////////////////////////////////////////////////
#ifdef WIN32
// there is a bug in the Microsoft VisualC++ compiler
// this class must be compiled with optimization off on Windows
# pragma optimize( "", off )
#endif
#include <TMCParticle.h>
#include <TFile.h>
#include <TSystem.h>
#include <TRandom.h>
#include <TROOT.h>
#include <TMath.h>
#include <TH2.h>
#include <TH3.h>
#include "ATLFast.h"
#include "ATLFMCMaker.h"
#include "ATLFTrackMaker.h"
#include "ATLFTrack.h"
const Float_t kPI = TMath::Pi();
const Float_t k2PI = 2*kPI;
const Int_t kMAXTRA = 100;
const Int_t kRECONS = BIT(16);
ClassImp(ATLFTrackMaker)
//_____________________________________________________________________________
ATLFTrackMaker::ATLFTrackMaker()
{
m_Ntracks = 0;
}
//_____________________________________________________________________________
ATLFTrackMaker::ATLFTrackMaker(const char *name, const char *title)
:ATLFMaker(name,title)
{
// Default Setters for tracks
SetMinPT();
SetMaxEta();
SetFieldType();
SetBLayer();
SetBeamConstraint();
m_Fruits = new TClonesArray("ATLFTrack",100, kFALSE);
m_BranchName = "Tracks";
m_Ntracks = 0;
// Please, how to do this optionally ??!!!
Save();
}
//_____________________________________________________________________________
ATLFTrackMaker::~ATLFTrackMaker()
{
//dummy
}
//_____________________________________________________________________________
ATLFTrack *ATLFTrackMaker::AddTrack(Int_t code, Int_t mcparticle)
{
// Add a new track to the list of tracks
//Note the use of the "new with placement" to create a new track object.
//This complex "new" works in the following way:
// tracks[i] is the value of the pointer for track number i in the TClonesArray
// if it is zero, then a new track must be generated. This typically
// will happen only at the first events
// If it is not zero, then the already existing object is overwritten
// by the new track parameters.
// This technique should save a huge amount of time otherwise spent
// in the operators new and delete.
TClonesArray &tracks = *(TClonesArray*)m_Fruits;
return new(tracks[m_Ntracks++]) ATLFTrack(code,mcparticle);
}
//_____________________________________________________________________________
void ATLFTrackMaker::Clear(Option_t *option)
{
// Reset Track Maker
m_Ntracks = 0;
ATLFMaker::Clear(option);
}
//_____________________________________________________________________________
void ATLFTrackMaker::Draw(Option_t *)
{
// Dummy Draw
}
//_____________________________________________________________________________
void ATLFTrackMaker::Init()
{
Int_t i,j,k;
// Tracks histograms
m_Mult = new TH1F("TraMult","tracks multiplicity",100,0,100);
// Smear Tracks
if (gATLFast->Smearing() == 1) {
// Open root file containing parameterised resolutions
//Text_t TrackResFileName[256];
//strcpy(TrackResFileName,gSystem->
// ExpandPathName("$ATLFAST/data/IDResolMu011.root"));
//tobe portable between Unix and NT, the following is required
Text_t *TrackResFileName = gSystem->ConcatFileName(gSystem->Getenv("ATLFAST"),"/data/IDResolMu011.root");
if ( gSystem->AccessPathName(TrackResFileName,kReadPermission) )
printf("Error: Resolution file not foundn");
// Change Directory
TDirectory *dirsav = gDirectory;
dirsav->ls();
TFile *file = new TFile(TrackResFileName,"READ");
if (TrackResFileName) delete [] TrackResFileName;
gROOT->cd();
// Bin Sizes
TH2F *TvIP = (TH2F*)file->Get("TvIP");
m_nEta = TvIP->GetXaxis()->GetNbins();
m_nPt = TvIP->GetYaxis()->GetNbins();
Float_t *EtaVal = new Float_t [m_nEta+2];
for (i=0; i<m_nEta+1; i++) {
EtaVal[i] = TvIP->GetXaxis()->GetBinLowEdge(i+1);
}
Float_t *PtVal = new Float_t [m_nPt+2];
for (i=0; i<m_nPt+1; i++) {
PtVal[i] = TvIP->GetYaxis()->GetBinLowEdge(i+1);
}
Float_t ParamN[10];
for (i=0; i<=9; i++) {
ParamN[i] = i;
}
// 3d Resolution Histograms
m_OP1 = new TH3F("m_OP1","Resolutions",m_nEta,EtaVal,m_nPt,PtVal,9,ParamN);
m_RP1 = new TH3F("m_RP1"," ",m_nEta,EtaVal,m_nPt,PtVal,9,ParamN);
// Read in Histograms
Char_t *ResNames[9] = {"TvIP","LtIP","Phi0","CoTanT","QdPT",
"InvPtPhi","InvPtA0","PhiA0","CoTTZ0"};
for (k=1; k<=9; k++) {
TH2F *H2 = (TH2F*)file->Get(ResNames[k-1]);
for (i=1; i<=m_nEta; i++) {
for (j=1; j<=m_nPt; j++) {
m_OP1->SetBinContent(m_OP1->GetBin(i,j,k),H2->GetCellContent(i,j));
}
}
}
delete [] EtaVal;
delete [] PtVal;
// Calculate resolution parameters
for (i=1; i<m_nEta; i++) {
for (j=1; j<m_nPt; j++) {
for (k=1; k<=5; k++) {
m_RP1->SetBinContent(m_RP1->GetBin(i,j,k),
TMath::Power(m_OP1->GetBinContent(m_OP1->GetBin(i,j,k)),2.0));
}
}
}
for (i=1; i<m_nEta; i++) {
for (j=1; j<m_nPt; j++) {
m_RP1->SetBinContent(m_RP1->GetBin(i,j,6),m_OP1->GetBinContent(m_OP1->GetBin(i,j,3))*
m_OP1->GetBinContent(m_OP1->GetBin(i,j,5))*
m_OP1->GetBinContent(m_OP1->GetBin(i,j,6)));
m_RP1->SetBinContent(m_RP1->GetBin(i,j,7),m_OP1->GetBinContent(m_OP1->GetBin(i,j,1))*
m_OP1->GetBinContent(m_OP1->GetBin(i,j,5))*
m_OP1->GetBinContent(m_OP1->GetBin(i,j,7)));
m_RP1->SetBinContent(m_RP1->GetBin(i,j,8),m_OP1->GetBinContent(m_OP1->GetBin(i,j,1))*
m_OP1->GetBinContent(m_OP1->GetBin(i,j,3))*
m_OP1->GetBinContent(m_OP1->GetBin(i,j,8)));
m_RP1->SetBinContent(m_RP1->GetBin(i,j,9),m_OP1->GetBinContent(m_OP1->GetBin(i,j,2))*
m_OP1->GetBinContent(m_OP1->GetBin(i,j,4))*
m_OP1->GetBinContent(m_OP1->GetBin(i,j,9)));
}
}
// Close file
delete file;
if (dirsav) dirsav->cd();
} // End gATLFast->Smearing(1) loop
}
//_____________________________________________________________________________
void ATLFTrackMaker::Finish()
{
// Function called by ATLFast::Finish at the end of the job
}
//_____________________________________________________________________________
Int_t ATLFTrackMaker::Make()
{
//.............................................
// This function searches for tracks, by scanning through
// the list of MC particles. If a track is found, its momentum is
// smeared using function RESMUO.
// Three options for track-momentum smearing are available:
// stand-alone Track System
// Inner Detector alone
// combined
// The parametrization for the momentum smearing is coded in RESMUO.
// Isolated and non-isolated tracks are stored in The TClonesArray of tracks
// and the energy clusters associated with them are removed.
// Tracks outside the ETA-coverage or below the p_T-threshold are lost.
//.............................................
Int_t i, k, km, KS, KF;
//........................................................
//... Get pointers to MCparticles arrays and ClonesArray
//........................................................
ATLFMCMaker *mcarlo = gATLFast->MCMaker();
TClonesArray *particles = mcarlo->Fruits();
TMCParticle *part, *parent;
Float_t *pPT = mcarlo->PT();
Float_t *pETA = mcarlo->Eta();
Int_t *Index= mcarlo->Index();
//.............................................
//...Loop over all particles. Find charged tracks
// traS: holds smeared track parameters
//.............................................
ATLFTrack *track;
Int_t trackMother[kMAXTRA], trackKF[kMAXTRA];
Float_t vert[3], pvert[4], tra[5], traS[5];
Float_t charge;
Float_t a0, z0, phi, cot, ptinv;
a0 = z0= phi= cot= ptinv= 0;
Float_t a0Crude, z0Crude, phiCrude, cotCrude, ptinvCrude;
Float_t corr[25];
m_Ntracks = 0;
if (!gATLFast->TrackFinding()) return 0;
//.............................................
//...for charged tracks calculate smeared and unsmeared parameters
//...store mother information
//.............................................
for (Int_t ind=0;ind<mcarlo->Nstables();ind++) {
i = Index[ind];
part = (TMCParticle*)particles->UncheckedAt(i);
KS = part->GetKS();
KF = part->GetKF();
if (KS <= 0 || KS > 10) continue;
if (pPT[i] <= m_MinPT) continue;
if (fabs(pETA[i]) > m_MaxEta) continue;
if (mcarlo->Charge(KF) == 0) continue;
trackKF[0] = KF;
trackMother[0] = part->GetParent();
Int_t nmother=0;
for (k=1;k<=6;k++) {
km = trackMother[k-1];
parent = (TMCParticle*)particles->UncheckedAt(km-1);
if (!parent) break;
nmother++;
trackMother[k] = parent->GetParent();
trackKF[k] = parent->GetKF();
if (trackMother[k] == 92) break;
}
charge = mcarlo->Charge(KF)/3;
//.............................................
// V is in mm, and we need it in cm, so divide by 10
//.............................................
vert[0] = 0.1*part->GetVx();
vert[1] = 0.1*part->GetVy();
vert[2] = 0.1*part->GetVz();
pvert[0] = part->GetPx();
pvert[1] = part->GetPy();
pvert[2] = part->GetPz();
pvert[3] = part->GetEnergy();
//.............................................
//.... reconstruct tracks
//.............................................
HelixParameters(charge, vert, pvert, tra);
a0Crude = tra[0];
z0Crude = tra[1];
phiCrude = tra[2];
cotCrude = tra[3];
ptinvCrude = tra[4];
//.............................................
//.... smear reconstructed tracks
//.............................................
if (gATLFast->Smearing() == 1) {
// Sigma array
TMatrix Sigma(5,5);
// Check for particle type and pt/eta bounds
if ( KF == 11 ) { // No Electron smearing yet
a0 = tra[0];
z0 = tra[1];
phi = tra[2];
cot = tra[3];
ptinv = tra[4];
} else
if ( KF == 13 ) { // Muon....gaussian only)
//Calcutate Sigmas
SmearParameters(pETA[i],pPT[i],tra,traS,Sigma);
a0 = traS[0];
z0 = traS[1];
phi = traS[2];
cot = traS[3];
ptinv = traS[4];
} else { // All others as Pions
//Calcutate Sigmas
SmearParameters(pETA[i],pPT[i],tra,traS,Sigma);
a0 = traS[0];
z0 = traS[1];
phi = traS[2];
cot = traS[3];
ptinv = traS[4];
}
} else if(gATLFast->Smearing() == 0) {
a0 = tra[0];
z0 = tra[1];
phi = tra[2];
cot = tra[3];
ptinv = tra[4];
}
//........................................................
//....Store track in Root ClonesArray
//........................................................
track = AddTrack(KF,i);
/// this is not very elegant, sorry
if(nmother > 0)track->SetMother(0,trackKF[0],trackMother[0]);
if(nmother > 1)track->SetMother(1,trackKF[1],trackMother[1]);
if(nmother > 2)track->SetMother(2,trackKF[2],trackMother[2]);
if(nmother > 3)track->SetMother(3,trackKF[3],trackMother[3]);
if(nmother > 4)track->SetMother(4,trackKF[4],trackMother[4]);
if(nmother > 5)track->SetMother(5,trackKF[5],trackMother[5]);
track->SetCrudeTrack(a0Crude, z0Crude, phiCrude, cotCrude, ptinvCrude);
track->SetSmearedTrack(a0, z0, phi, cot, ptinv);
track->SetCorrelations(corr);
// mark original particle as being reconstructed
part->SetBit(kRECONS);
}
m_Mult->Fill(m_Ntracks+0.01);
return 0;
}
//_____________________________________________________________________________
void ATLFTrackMaker::PrintInfo()
{
ATLFMaker::PrintInfo();
}
//_____________________________________________________________________________
void ATLFTrackMaker::HelixParameters(Float_t charge, Float_t *vert1, Float_t *pvert1, Float_t *b)
{
// Returns helix parameters in polar coordinates from the original position
// and momentum at the vertex of the track.
// Input Arguments:
// ================
//
// charge = particle charge
// vert = vertex coordinates
// pvert = quadrimpulse of the particle
//
// Output Arguments:
// =================
//
// b[0] Impact parameter ! output is ATLAS standards
// b[1] Z of perigee ! in cm
// b[2] phi of helix (-pi<phi<pi) ! in rad
// b[3] cot(theta)
// b[4] (1/pt)*charge ! in 1/GeV
//
// Author:
// =======
//
// Tarta
//======================================================================
Int_t i;
Float_t q, bmagnt, dotsign;
Float_t c, d, phi0, cote, z0,x0,y0,x01,y01;
Double_t dd,cc,pphi0,zz0,rrho,ccote,xx0,yy0,px,py,pz,dot;
Double_t xv,yv,zv,cosa,sina,pt,diffv,alf,bet,angdif;
Double_t vert[3], pvert[4];
for (i=0;i<3;i++) vert[i] = vert1[i];
for (i=0;i<4;i++) pvert[i] = pvert1[i];
bmagnt = 2;
px = pvert[0];
py = pvert[1];
pz = pvert[2];
q = charge;
pt = TMath::Sqrt(px*px+py*py);
xv = vert[0];
yv = vert[1];
zv = vert[2];
// Get track polar coordinate parameters
cc = q*0.00149898*bmagnt/pt;
c = cc;
rrho = q/(2*cc);
sina = py/pt;
cosa = px/pt;
xx0 = xv + q*rrho*sina;
yy0 = yv - q*rrho*cosa;
x0 = xx0;
y0 = yy0;
dd = q*(TMath::Sqrt(xx0*xx0+yy0*yy0) - rrho);
d = -dd;
alf = TMath::ATan2(yy0,xx0);
if (TMath::Abs(yv) < 1e-4 && TMath::Abs(xv) < 1e-4) {
pphi0 = TMath::ATan2(py,px);
} else {
bet = TMath::ATan2(yv,xv);
angdif = bet - alf;
if (angdif > kPI) angdif -= k2PI;
else if(angdif < -kPI) angdif += k2PI;
if (angdif > 0) pphi0 = alf + 0.5*kPI;
else pphi0 = alf - 0.5*kPI;
}
x01 = TMath::Sin(pphi0)*(1 + 2*c*dd)/(2*c);
y01 = -TMath::Cos(pphi0)*(1 + 2*c*dd)/(2*c);
if (x01*x0 < 0 || y01*y0 < 0) pphi0 += kPI;
if (pphi0 < 0) pphi0 += k2PI;
else if(pphi0 > k2PI) pphi0 -= k2PI;
phi0 = pphi0;
ccote = pz/pt;
cote = ccote;
// Protection in case of very small impact parameters
diffv = xv*xv +yv*yv -dd*dd;
if (diffv < 0) diffv = 0;
dot = xv*px + yv*py;
if (dot < 0) dotsign = -1;
else dotsign = 1;
zz0 = zv -dotsign*ccote*TMath::ASin(cc*TMath::Sqrt(diffv/(1+2*cc*dd)))/cc;
z0 = zz0;
// output
b[0] = d;
b[1] = z0;
b[2] = phi0;
b[3] = cote;
b[4] = q/pt ;
}
//_____________________________________________________________________________
void ATLFTrackMaker::Resolution(Float_t eta, Float_t pt, Double_t *Sigsq)
{
Int_t i;
eta = TMath::Abs(eta);
Int_t bineta = m_RP1->GetXaxis()->FindBin(eta);
Int_t binpt = m_RP1->GetYaxis()->FindBin(pt);
for (i=1; i<=9; i++) {
// Interpolate between bineta+1,binpt+1 and bineta,binpt
Float_t res11 = m_RP1->GetBinContent(m_RP1->GetBin(bineta,binpt,i));
Float_t res21 = m_RP1->GetBinContent(m_RP1->GetBin(bineta+1,binpt,i));
Float_t res12 = m_RP1->GetBinContent(m_RP1->GetBin(bineta,binpt+1,i));
Float_t res22 = m_RP1->GetBinContent(m_RP1->GetBin(bineta+1,binpt+1,i));
Float_t Li1 = res11 + ( res21 - res11 )*
( (eta - (m_RP1->GetXaxis()->GetBinLowEdge(bineta))) /
((m_RP1->GetXaxis()->GetBinLowEdge(bineta+1)) -
(m_RP1->GetXaxis()->GetBinLowEdge(bineta))));
Float_t Li2 = res12 + ( res22 - res12 )*
( (eta - (m_RP1->GetXaxis()->GetBinLowEdge(bineta))) /
((m_RP1->GetXaxis()->GetBinLowEdge(bineta+1)) -
(m_RP1->GetXaxis()->GetBinLowEdge(bineta))));
// Interpolate between Li1 and Li1
Float_t asq = (Li2*TMath::Power(m_RP1->GetYaxis()->GetBinLowEdge(binpt+1),2) -
Li1*TMath::Power(m_RP1->GetYaxis()->GetBinLowEdge(binpt),2)) /
(TMath::Power(m_RP1->GetYaxis()->GetBinLowEdge(binpt+1),2) -
TMath::Power(m_RP1->GetYaxis()->GetBinLowEdge(binpt),2));
Float_t bsq = (Li2 - Li1)*(TMath::Power(m_RP1->GetYaxis()->GetBinLowEdge(binpt)*
m_RP1->GetYaxis()->GetBinLowEdge(binpt+1),2))/
(TMath::Power(m_RP1->GetYaxis()->GetBinLowEdge(binpt),2) -
TMath::Power(m_RP1->GetYaxis()->GetBinLowEdge(binpt+1),2));
Float_t sigsq = asq + bsq/TMath::Power(pt,2);
Sigsq[i-1] = sigsq;
}
}
//_____________________________________________________________________________
void ATLFTrackMaker::SmearParameters(Float_t eta, Float_t pt, Float_t*tra,
Float_t *traS, TMatrix &Sigma)
{
// Calculates Smeared parameters for a track
//
// Input Arguments:
// ================
//
// Eta and Pt of track to be smeared.
// Unsmeared Track parameters
// Output Arguments:
// =================
// Smeared track parameters
// Sigma correlation matrix
//
//======================================================================
Int_t i;
Float_t conv = 10000;
Float_t convR = 1000;
Double_t X[5] = {0,0,0,0,0}; // Vector of correlated gaussian variables
Double_t SigSq[9] = {0,0,0,0,0,0,0,0,0};
// Fill SigSq 1D Histo with Sigma Parameters
Resolution(eta,pt,SigSq);
// Fill Sigma Matrix
// Diagonal elements
for (i = 0; i < 5; i++){
Sigma(i,i) = SigSq[i];
}
// Covariance terms
Sigma(0,2) = SigSq[7];
Sigma(0,4) = SigSq[6];
Sigma(1,3) = SigSq[8];
Sigma(2,4) = SigSq[5];
Sigma(2,0) = Sigma(0,2);
Sigma(4,0) = Sigma(0,4);
Sigma(3,1) = Sigma(1,3);
Sigma(4,2) = Sigma(2,4);
// ==> call routines to generate correlated gaussian random numbers
TMatrix Cigma(5,5);
dcorset(Sigma,Cigma,5);
dcorgen(Cigma,X,5);
for (i=0; i<2; i++ )
traS[i] = (tra[i]*conv + X[i])/conv;
traS[2] = (tra[2]*convR + X[2])/convR;
if (traS[2] > k2PI )
traS[2] -= k2PI;
else
if (traS[2] < 0 )
traS[2] += k2PI;
for(i=3; i<5; i++)
traS[i] = tra[i] + X[i]/1000;
}
//______________________________________________________________________________
void ATLFTrackMaker::dcorset(TMatrix &sigma, TMatrix &cigma, Int_t n)
{
Int_t i,j,k;
Double_t sum;
// Compute square root of matrix sigma
for (i=0; i<n; i++) {
sum = 0;
for (j=0; j<i; j++) {
sum += cigma(i,j)*cigma(i,j);
}
cigma(i,i) = sqrt(TMath::Abs(sigma(i,i) - sum));
// Off Diagonal terms
for (k=i+1; k<n; k++) {
sum = 0;
for (j=0; j<i; j++) {
sum += cigma(k,j)*cigma(i,j);
}
cigma(k,i) = (sigma(k,i) - sum)/cigma(i,i);
}
}
}
//______________________________________________________________________________
void ATLFTrackMaker::dcorgen(TMatrix &cigma, Double_t *x, Int_t n)
{
Int_t i,j;
Int_t nmax = 100;
if (n > nmax ) {
printf("Error in dcorgen: array overflown");
}
Float_t tmp;
for (i=0; i<n; i++) {
x[i] = 0.0;
for (j=0; j<=i; j++) {
tmp = gRandom->Gaus(0,1);
x[i] += cigma(i,j)*tmp;
}
}
}
//______________________________________________________________________________
void ATLFTrackMaker::Streamer(TBuffer &R__b)
{
// Stream an object of class ATLFTrackMaker.
if (R__b.IsReading()) {
Version_t R__v = R__b.ReadVersion();
ATLFMaker::Streamer(R__b);
R__b >> m_Ntracks;
R__b >> m_MinPT;
R__b >> m_MaxEta;
R__b >> m_Mult;
if ( R__v < 2 ) return;
R__b >> m_FieldType;
R__b >> m_BLayer;
R__b >> m_BeamConstraint;
} else {
R__b.WriteVersion(ATLFTrackMaker::IsA());
ATLFMaker::Streamer(R__b);
R__b << m_Ntracks;
R__b << m_MinPT;
R__b << m_MaxEta;
R__b << (TObject*)m_Mult;
R__b << m_FieldType;
R__b << m_BLayer;
R__b << m_BeamConstraint;
}
}
//______________________________________________________________________________
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.