/* * PndLmdCalStrip.CXX * * Created on: Jul 28, 2009 * Author: huagen */ //******************************************************************** //the implementation of PndLmdCalStrip class //******************************************************************** #include "PndLmdCalStrip.h" #include "TMath.h" #include "TRandom3.h" #include #include #include "PndLmdStrips.h" #include "PndLmdDigiPara.h" using namespace std; ClassImp(PndLmdCalStrip); PndLmdCalStrip::PndLmdCalStrip() { fPitch = 0; fOrient = 0; fAnchor = TVector2(0,0); fThreshold = 0; fNoise = 0; fNrFEChannels = 0; fNrStrips = 0; } PndLmdCalStrip::PndLmdCalStrip(Double_t pitch,Double_t orient, Int_t nrStrips, Int_t nrFeChannels, const TVector2& firstStripAnchor,Double_t threshold, Double_t noise,Double_t sigma) //:fPitch(pitch),fOrient(orient),fNrStrips(nrStrips),fNrFEChannels(nrFeChannels), // fAnchor(firstStripAnchor),fThreshold(threshold),fNoise(noise) { fPitch = pitch; fOrient = orient; fAnchor = firstStripAnchor; fThreshold = threshold; fNoise = noise; fNrFEChannels = nrFeChannels; fNrStrips = nrStrips; fSigma = sigma; fRNG = new TRandom3(); fStripDir.Set(cos(fOrient),sin(fOrient)); fOrthoDir.Set(sin(fOrient),cos(fOrient)); } PndLmdCalStrip::PndLmdCalStrip(const PndLmdDigiPara* digipara, SensorSegment segment) { if(segment == kLEFT) { fPitch = digipara->GetLeftPitch(); fOrient = digipara->GetLeftOrient(); fAnchor = digipara->GetLeftAnchor(); fNrStrips = digipara->GetNrLeftFe()*digipara->GetNrFeChannels(); fStripDir.Set(cos(fOrient),sin(fOrient)); fOrthoDir.Set(-sin(fOrient),cos(fOrient)); } else if(segment==kRIGHT) { fPitch = digipara->GetRightPitch(); fOrient = digipara->GetRightOrient(); fAnchor = digipara->GetRightAnchor(); fNrStrips = digipara->GetNrRightFe()*digipara->GetNrFeChannels(); fStripDir.Set(cos(fOrient),sin(fOrient)); fOrthoDir.Set(sin(fOrient),-cos(fOrient)); } fNrFEChannels = digipara->GetNrFeChannels(); fThreshold = digipara->GetThreshold(); fNoise = digipara->GetNoise(); fSigma = digipara->GetGausSigma(); fRNG = new TRandom3(); } /* //calculate the strips for curved strips PndLmdCalStrip::PndLmdCalStrip(const PndLmdDigiPara* digipara, SensorSide side) { if(side == kTOP) fPitch = digipara->GetCirclePitch(); // fOrient = digipara->GetRightOrient(); fAnchor = digipara->GetCircleAnchor(); fNrStrips = digipara->GetNrCircleFe()*digipara->GetNrFeChannels(); // fStripDir.Set(con(fOrient),sin(fOrient)); // fOrthoDir.Set(sin(fOrient),-cos(fOrient)); fNrFEChannels = digipara->GetNrFeChannels(); fThreshold = digipara->GetThreshold(); fNoise = digipara->GetNoise(); } */ std::vector PndLmdCalStrip::GetStrips(Double_t inx, Double_t iny, Double_t inz, Double_t outx, Double_t outy, Double_t outz,Double_t eLoss) { //2d projection of trajectory TVector2 in(inx,iny); TVector2 out(outx, outy); TVector2 tpath= out-in; Double_t path = tpath.Mod(); std::vector strips; // Double_t fract_Q; if(path<1e-18) { std::cout<<"-w- PndLmdCalStrip::GetStrips: No Trajectory inside sensor!"<Double_t(fNrStrips)&&nuOut>Double_t(fNrStrips)){ std::cout<<"-i- warning: Hit outside active area."<(Double_t(fNrStrips))){ Q *=((Double_t)fNrStrips-nuIn)/(nuOut-nuIn); nuIn = Double_t(fNrStrips); } //is the Out-point inside the active area? if(nuOut<0.){ Q *=(nuIn)/(-nuOut+nuIn); nuOut = 0; }else if(nuOut>(Double_t(fNrStrips))){ Q *=((Double_t)fNrStrips-nuIn)/(nuOut-nuIn); nuOut = Double_t(fNrStrips); } //only 1 hit happened if(inStrip==outStrip){ // fract_Q = Q; strips.push_back(PndLmdStrips(inStrip,path,Q)); // cout<<"Only one hit happened and the strip fired is "<fThreshold) strips.push_back(PndLmdStrips(inStrip,path1,Q1)); nrHits++; Q -= Q1; //calculate the portion of track in the last strip Int_t prevOut = Int_t(nuOut+0.5+0.5*dir); Double_t path2 = std::fabs(nuOut-prevOut)*dPath; Double_t Q2 = dQ * path2; if(Q2 > fThreshold) strips.push_back(PndLmdStrips(outStrip,path2,Q2)); nrHits++; Q -= Q2; //distribute the charge among the intermediate strips nextIn = Int_t(nextIn-0.5+0.5*dir); prevOut = Int_t(prevOut-0.5+0.5*dir); for(Int_t n = nextIn; n != prevOut; n +=Int_t(dir)){ if(dQ > fThreshold) strips.push_back(PndLmdStrips(n,dPath,dQ)); nrHits++; Q -= dQ; /* cout<<"more than 1 hits: the strips fired are :" <<"inStrip is"<1 loop return strips; } //the public method for strip digitization, which the charge diffusion was taken into account std::vector PndLmdCalStrip::DigiStripFromPoint(Double_t inx, Double_t iny, Double_t inz, Double_t outx, Double_t outy,Double_t outz,Double_t eLoss) { //call the GetStrips to initialize the strips std::vector strips = GetStrips(inx, iny, inz, outx, outy,outz,eLoss); cout<<" Number of strips initially fired :"<=fThreshold) strips.push_back(PndLmdStrips(index,path_real,dQ)); // cout<<"the dQ for the center strip is "<=fThreshold){ //charge distribution at the head of the track charge1 = ChargeDiffusion((i-1)*fPitch,Path1,path,dir,Q);//charge distribution from -infinite to strip of index-1 charge2 = ChargeDiffusion(i*fPitch,Path1,path,dir,Q); q_head = charge2 - charge1; id = index + i-1; dQ = AddNoise(q_head); if(dQ>fThreshold) strips.push_back(PndLmdStrips(id,0,dQ)); } if(dQR>=fThreshold){ //charge distribution at the head of the track charge1 = ChargeDiffusion(j*fPitch,Path1,path,dir,Q);//charge distribution from -infinite to strip of index-1 charge2 = ChargeDiffusion((j+1)*fPitch,Path1,path,dir,Q); q_tail = charge2 - charge1; id = index + j; dQ = AddNoise(q_tail); if(dQ>fThreshold) strips.push_back(PndLmdStrips(id,0,dQ)); } //the charge left on the center strip Double_t restQ = ChargeDiffusion(fPitch,Path1,path,dir,Q) - ChargeDiffusion(0,Path1,path,dir,Q); strips.push_back(PndLmdStrips(index,path_real,AddNoise(restQ))); } }//for size == 1 // else for the size>1, means the number of strips fired greater than 1 else{ //for the strips Double_t head_Path, tail_Path; std::vector strip_head; std::vector strip_tail; std::vector strip_middle; PndLmdStrips str_head; PndLmdStrips str_tail; PndLmdStrips str_middle; if(inStrip0;i--) {str_middle = strips.at(i);} } strips.erase(strips.begin(),strips.begin()+size); //strips initialization Int_t index_head = str_head.GetIndex(); Double_t path_head = str_head.GetPath(); //Double_t q_head = str_head.GetCharge(); Int_t index_tail = str_tail.GetIndex(); Double_t path_tail = str_tail.GetPath(); // Double_t q_tail = str_tail.GetCharge(); //for the left strip of the head strip Double_t charge1, charge2, dQ; Int_t id; Double_t ipath; //Double_t q_head,q_tail; Double_t dq_head = ChargeDiffusion(0,head_Path,path,dir,total_Q);//the diffusion charge from -infinite to index Double_t dq_tail = total_Q- ChargeDiffusion(size*fPitch,tail_Path,path,(-1)*dir,total_Q);//the diffusion charge from index+size to infinite Double_t q_head = ChargeDiffusion(fPitch,head_Path,path,dir,total_Q);// the charge left on the strip index Double_t q_tail = ChargeDiffusion(size*fPitch,tail_Path,path,(-1)*dir,total_Q) -ChargeDiffusion((size-1)*fPitch,tail_Path,path,(-1)*dir,total_Q); Double_t dQ_head = AddNoise(dq_head); Double_t dQ_tail = AddNoise(dq_tail); Double_t Q_head = AddNoise(q_head); Double_t Q_tail = AddNoise(q_tail); if(dQ_head>fThreshold) { id = index_head-1; strip_head.push_back(PndLmdStrips(id,0,dQ_head));} if(Q_head>fThreshold) { id = index_head; strip_head.push_back(PndLmdStrips(id,path_head,Q_head));} strips.insert(strips.begin(),strip_head.begin(),strip_head.end()); if(dQ_tail>fThreshold) { id = index_tail+1; strip_tail.push_back(PndLmdStrips(id,0,dQ_tail));} if(Q_tail>fThreshold) { id = index_tail; strip_head.push_back(PndLmdStrips(id,path_tail,Q_tail));} strips.insert(strips.end(),strip_tail.begin(),strip_tail.end()); //for the strips between the head and the tail for(Double_t i= 1;i< size-1;i++){ charge1 = ChargeDiffusion(i*fPitch,head_Path,path,dir,total_Q); charge2 = ChargeDiffusion((i+1)*fPitch,head_Path,path,dir,total_Q); dQ = charge2 - charge1; std::vector::iterator iter_ii=strip_middle.begin(); id = iter_ii->GetIndex(); ipath = iter_ii->GetPath(); if(dQ>fThreshold) strip_middle.push_back(PndLmdStrips(id,ipath,dQ)); iter_ii++; }// strips.insert(strips.begin()+strip_head.size(),strip_middle.begin(),strip_middle.begin()+strip_middle.size()); }//for size>1 return strips; }//digitize the strips Double_t PndLmdCalStrip::ChargeFromEloss(Double_t eloss) const {return eloss/3.61*1e9;} //calculator the strips from the anchor point Double_t PndLmdCalStrip::CalStripFromPoint(Double_t x, Double_t y) { return ((x-fAnchor.X())*fOrthoDir.X()+(y-fAnchor.Y())*fOrthoDir.Y())/fPitch;} Double_t PndLmdCalStrip::CalCurveStripFromPoint(Double_t x, Double_t y) {return (std::sqrt((x-fAnchor.X())*(x-fAnchor.X())+(y-fAnchor.Y())*(y-fAnchor.Y()))-30.0)/fPitch;} Int_t PndLmdCalStrip::CalChannelFromStrip(Int_t stripNr) {return stripNr%fNrFEChannels;} Int_t PndLmdCalStrip::CalFeFromStrip(Int_t stripNr) {return stripNr/fNrFEChannels;} //add the noise to the signal with gaussian distribution Double_t PndLmdCalStrip::AddNoise(Double_t charge) { Double_t smeared = fRNG->Gaus(charge,fNoise); return smeared; } //charge diffusion during collecting after ionization, //here, the u means the range which approaching the mean value a. //the area was expressed by the Erf function as below // f(x) = 1/2*(1+TMath::Erf((u-a)/sigma*TMath::sqrt(2.)))) Double_t PndLmdCalStrip::ChargeDiffusion(Double_t u,Double_t a,Double_t path,Double_t dir,Double_t Q) { Int_t N = 150; Double_t Npath = path/N; Double_t NQ = Q/N; Double_t Delta; //the area of Gaussian from (-infinite,u) Double_t dQ=0; //the Q in the given range if(path==0.) { Delta = 0.5 *(1+TMath::Erf((u-a)/fSigma*std::sqrt(2.))); dQ = NQ*Delta; } else { if(dir>0.) { for(Double_t da=a;da<=a+path;da+=Npath) { Delta = 0.5*(1+TMath::Erf((u-da)/fSigma*std::sqrt(2.))); dQ +=NQ * Delta; } } if(dir<0) { for(Double_t da=a-path;da<=a;da+=Npath) { Delta = 0.5*(1+TMath::Erf((u-da)/fSigma*std::sqrt(2.))); dQ +=NQ * Delta; } }//dir<<0 }//path>0 return dQ; } void PndLmdCalStrip::print() const { std::cout<<"-I- PndLmdCalStrip Info :"<