Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

DedxCalibDipAngle Class Reference

#include <DedxCalibDipAngle.h>

Inheritance diagram for DedxCalibDipAngle:

DedxCalib DedxCalib List of all members.

Public Member Functions

void AnalyseHists ()
void AnalyseHists ()
void BookHists ()
void BookHists ()
void checkSelections ()
void checkSelections ()
 DedxCalibDipAngle (const std::string &name, ISvcLocator *pSvcLocator)
 DedxCalibDipAngle (const std::string &name, ISvcLocator *pSvcLocator)
virtual StatusCode execute ()
virtual StatusCode execute ()
void FillHists ()
void FillHists ()
void FillHistsFromTree ()
void FillHistsFromTree ()
void FillTestHists ()
void FillTestHists ()
void FillTree ()
void FillTree ()
virtual StatusCode finalize ()
virtual StatusCode finalize ()
virtual Float_t GetChargeOffCorr (Int_t il, Int_t iw, Float_t ChargeIn, Float_t path, Float_t doca, Float_t EAngle, Float_t dipAngle, Float_t z)
virtual Float_t GetChargeOffCorr (Int_t il, Int_t iw, Float_t ChargeIn, Float_t path, Float_t doca, Float_t EAngle, Float_t dipAngle, Float_t z)
Float_t GetChargeOnCorr (Int_t il, Int_t iw, Float_t ChargeIn, Float_t path, Float_t doca, Float_t EAngle, Float_t dipAngle, Float_t z)
Float_t GetChargeOnCorr (Int_t il, Int_t iw, Float_t ChargeIn, Float_t path, Float_t doca, Float_t EAngle, Float_t dipAngle, Float_t z)
virtual StatusCode initialize ()
virtual StatusCode initialize ()
void ReadBetheBlochParameters ()
void ReadBetheBlochParameters ()
void ReadCalibdEdxParameters ()
void ReadCalibdEdxParameters ()
void WriteCalibdEdxParameters ()
void WriteCalibdEdxParameters ()
void WriteHists ()
void WriteHists ()

Public Attributes

double m_par [3]

Protected Attributes

int calib_flag
bool ddgflag
IDedxCorrecSvc * exsvc
IDedxCorrecSvc * exsvc
IMdcGeomSvc * geosvc
IMdcGeomSvc * geosvc
bool ggsflag
bool layergflag
MsgStream log
std::string m_constrootfile
int m_correc_flag
std::string m_inputfile
int m_par_flag
int m_phshape_flag
std::string m_rootfile
bool wiregflag
bool zdepflag

Private Attributes

TObjArray * m_caliblist
TObjArray * m_caliblist
TRandom * m_rand
TRandom * m_rand
TObjArray * m_satlist
TObjArray * m_satlist

Friends

void fcnQNew (Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
void fcnQNew (Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
void fcnQOLD (Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
void fcnQOLD (Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)

Constructor & Destructor Documentation

DedxCalibDipAngle::DedxCalibDipAngle (  const std::string &  name,
ISvcLocator *  pSvcLocator
) 
  

00056   : DedxCalib(name, pSvcLocator)
00057 {
00058   log<<MSG::INFO<<"DedxCalibDipAngle::DedxCalibDipAngle( )..."<<endreq;
00059 }

DedxCalibDipAngle::DedxCalibDipAngle (  const std::string &  name,
ISvcLocator *  pSvcLocator
) 
 

Member Function Documentation

void DedxCalibDipAngle::AnalyseHists (   )  [virtual]
 

Implements DedxCalib.

void DedxCalibDipAngle::AnalyseHists (   )  [virtual]
 

Implements DedxCalib.

00244 { 
00245   log<<MSG::INFO<<"DedxCalibDipAngle::AnalyseHists( )..."<<endreq;
00246    int long num_Qraw =0;
00247    double averdE_Fill;
00248    std::vector<Float_t>* vde_dedx=0;
00249    for(Int_t i=0;i<NUMTrk;i++)
00250     {
00251       //int P_Bin=(Int_t)floor((LOGP[i]-BetaGammaLogMinSat)/stepP);
00252       int P_Bin=(Int_t)floor((LOGP[i]-MomentMinSat)/stepP);
00253       int SinDip_Bin =(Int_t)floor((fabs(SINDip[i])-DipAngMin)/DipAngStep);
00254       vde_dedx=new std::vector<Float_t> (NUMHits[i]);
00255       for(Int_t j=0;j<NUMHits[i];j++)
00256         {
00257           (*vde_dedx)[j]=QRaw[num_Qraw];
00258           num_Qraw++;
00259        }
00260       std::stable_sort(vde_dedx->begin(),vde_dedx->end());
00261       //int num_Samples=(Int_t)(0.8*NUMHits/[i]);
00262       int num_Samples=(Int_t)floor(0.8*NUMHits[i]);
00263       for(Int_t ns=0;ns<num_Samples;ns++) averdE_Fill+=(*vde_dedx)[ns];
00264       averdE_Fill/=num_Samples;
00265       QvsDipAng_P[P_Bin][SinDip_Bin] ->Fill(averdE_Fill); 
00266       averdE_Fill=0;
00267       delete vde_dedx;
00268     }
00269 
00270   Int_t numpar;
00271   if(OldAnalyse ==1) numpar=3;
00272   else numpar=8;  
00273   MinuitdE = new TMinuit(numpar);
00274   if( OldAnalyse ==1){
00275     MinuitdE->SetFCN(fcnQOLD);
00276   }
00277   else{
00278     MinuitdE->SetFCN(fcnQNew);
00279   }
00280 
00281   Double_t* arglist= new Double_t [numpar*2];
00282   Int_t    ierflg = 0;
00283   arglist[0]=-1;
00284   MinuitdE->SetMaxIterations(1000000);
00285   MinuitdE->mnexcm("SET PRINT", arglist,1,ierflg);
00286   arglist[0] = 0;
00287   MinuitdE->mnexcm("SET NOW", arglist,0,ierflg);
00288   static Double_t vstart[10]; 
00289   if(UseOldFunctions==1)
00290     {
00291       vstart[0]=1.0;
00292       vstart[1]=1.0;
00293       vstart[2]=1.0;
00294 
00295       //vstart[0]=CoupalParOld[0]*1e+04;
00296       //vstart[1]=CoupalParOld[1]*1e+01;
00297       //vstart[2]=CoupalParOld[2]*1e+09;
00298     }
00299   else
00300     {
00301       vstart[0]=1.0;
00302       vstart[1]=1.0;
00303       vstart[2]=1.0;
00304       //std::vector<Float_t>*   vecsat=CalibdEdxParameters->GetVecSat();
00305       //vstart[0]=(*vecsat)[0]*1e+04;
00306       //vstart[1]=(*vecsat)[1]*1e+01;
00307       //vstart[2]=(*vecsat)[2]*1e+09;
00308     }
00309   log<<MSG::INFO<<"\n parSatQ["<<0<<"]="<<setprecision(11)<<vstart[0]<<endreq;
00310   log<<MSG::INFO<<"\n parSatQ["<<1<<"]="<<setprecision(11)<<vstart[1]<<endreq;
00311   log<<MSG::INFO<<"\n parSatQ["<<2<<"]="<<setprecision(11)<<vstart[2]<<endreq;
00312   
00313   vstart[3]=1.0;
00314   vstart[4]=1.0;
00315   vstart[5]=1.0;
00316   vstart[6]=1.0;
00317   vstart[7]=1.0;
00318 
00319   //vstart[3]=DipAngPar[0];
00320   //vstart[4]=DipAngPar[1];
00321   //vstart[5]=DipAngPar[2];
00322   //vstart[6]=DipAngPar[3];
00323   //vstart[7]=DipAngPar[4];
00324 //    vstart[8]=DipAngPar[5];
00325   //    vstart[9]=DipAngPar[6];
00326  
00327   static Double_t step[10] = {0.1,0.05,0.01,0.01,0.005,0.005,0.005,0.005,0.005,0.005};
00328   //if(_FixDipAngPar.value()==true)
00329   //  {
00330   //    for(Int_t ia=3;ia<8;ia++)
00331   //      {
00332   //        step[ia]=0;
00333   //      }
00334   //  }
00335   MinuitdE->mnexcm("SET ERR", arglist,1,ierflg);
00336   MinuitdE->DefineParameter(0,"p0",vstart[0],step[0],0.001,10);
00337   MinuitdE->DefineParameter(1,"p1",vstart[1],step[1],0.001,10);
00338   //       MinuitdE->DefineParameter(2,"p2",vstart[2],step[2],0.0,30);
00339   MinuitdE->DefineParameter(2,"p2",0,step[2],0.0,30);
00340 
00341   if(OldAnalyse==0)
00342     {
00343       MinuitdE->DefineParameter(3,"p3",vstart[3],step[3],0.85,1.15);
00344       MinuitdE->DefineParameter(4,"p4",vstart[4],step[4],-0.1,0.1);
00345       MinuitdE->DefineParameter(5,"p5",vstart[5],step[5],-0.1,0.1);
00346       MinuitdE->DefineParameter(6,"p6",vstart[6],step[6],-0.1,0.1);
00347       MinuitdE->DefineParameter(7,"p7",vstart[7],step[7],-0.1,0.1);
00348       // MinuitdE->DefineParameter(8,"p8",vstart[8],step[8],-0.1,0.1);
00349       // MinuitdE->DefineParameter(9,"p9",vstart[9],step[9],-0.1,0.1);
00350     }
00351 
00352   MinuitdE->mnexcm("HESSE",  arglist,0,ierflg);
00353   //     MinuitdE->mnexcm("MINOs 1000 3 3",   arglist,0,ierflg);
00354   MinuitdE->mnexcm("MIGRAD",arglist,numpar,ierflg);
00355   MinuitdE->mnimpr();
00356   Double_t amin,edm,errdef;
00357   Int_t nvpar,nparx,icstat;
00358   Double_t*  par=new Double_t [numpar] ;
00359   Double_t*  parerr=new Double_t [numpar] ;
00360   MinuitdE->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
00361   MinuitdE->mnprin(numpar,amin);
00362 
00363   for(Int_t i=0;i<numpar;i++)
00364     {
00365       MinuitdE->GetParameter(i,par[i],parerr[i]);
00366     }
00367   if(OldAnalyse==0)
00368     {
00369       log<<MSG::INFO<< " p0: " <<setprecision(12) << par[0]*1e-04 << "\n "
00370                     << " p1: " << setprecision(12) <<par[1]*1e-01 << "\n "
00371                     << " p2: " << setprecision(12) <<par[2]*1e-09 << endreq;
00372  
00373       MinuitdE->DefineParameter(0,"p0",par[0],0,0.01,10);
00374       MinuitdE->DefineParameter(1,"p1",par[1],0,0.01,10);
00375       MinuitdE->DefineParameter(2,"p2",par[2],0.01,0.0,20);
00376       //             MinuitdE->DefineParameter(2,"p2",0,0.0,0.0,20);
00377  
00378       MinuitdE->DefineParameter(3,"p3",par[3],0,0.85,1.15);
00379       MinuitdE->DefineParameter(4,"p4",par[4],0,-0.1,0.1);
00380       MinuitdE->DefineParameter(5,"p5",par[5],0,-0.1,0.1);
00381       MinuitdE->DefineParameter(6,"p6",par[6],0,-0.1,0.1);
00382       MinuitdE->DefineParameter(7,"p7",par[7],0,-0.1,0.1);
00383       //                 MinuitdE->DefineParameter(8,"p6",par[8],0,-0.1,0.1);
00384       //                 MinuitdE->DefineParameter(9,"p7",par[9],0,-0.1,0.1);
00385  
00386       MinuitdE->mnexcm("MIGRAD",arglist,numpar,ierflg);
00387       MinuitdE->mnimpr();
00388       for(Int_t i=0;i<numpar;i++)
00389         {
00390           MinuitdE->GetParameter(i,par[i],parerr[i]);
00391         }
00392     }
00393   log<<MSG::INFO<< " p0: " <<setprecision(12) << par[0]*1e-04 << "\n "
00394                 << " p1: " << setprecision(12) <<par[1]*1e-01 << "\n "
00395                 << " p2: " << setprecision(12) <<par[2]*1e-09 << endreq;
00396   if(OldAnalyse==0)
00397     {
00398       for(Int_t i=3;i<numpar;i++)
00399         {
00400           log<<MSG::INFO << " a[ " <<i<<"]:"<<setprecision(12) << par[i] << endreq;
00401         }
00402     }
00403  /* std::vector<Float_t>* vs=CalibdEdxParameters->GetVecSat();
00404   (*vs)[0]=par[0]*1e-04;
00405   (*vs)[1]=par[1]*1e-01;
00406   (*vs)[2]=par[2]*1e-09;
00407   if(_SaveCalibdEdxParameters.value()==true)
00408     {
00409       CalibdEdxParameters->SetVecSat(*vs);
00410     }
00411   else
00412     {
00413       log<<MSG::INFO<<"\n result has not been saved !"<<endreq;
00414     }
00415 */
00416   if(OldAnalyse==1) {
00417      for(Int_t i=0; i<3; i++)
00418      {SavePar_old[i] = par[i];
00419       //cout<<"SavePar_old[ "<<i<<" ] ="<<SavePar_old[i]<<endl; 
00420      } 
00421   }
00422   else{
00423      for(Int_t i=0; i<8; i++)
00424      {SavePar_new[i] = par[i];}
00425   }
00426   
00427   delete [] arglist;
00428   delete [] par;
00429   delete [] parerr;
00430 
00431 
00432 
00433 }

void DedxCalibDipAngle::BookHists (   )  [virtual]
 

Implements DedxCalib.

void DedxCalibDipAngle::BookHists (   )  [virtual]
 

Implements DedxCalib.

00105 {
00106   log<<MSG::INFO<<"DedxCalibDipAngle::BookHists( )..."<<endreq;  
00107   char NameHist[40];
00108   QvsDipAngP=new  TProfile* [NumSlicesBGSat];
00109   QvsDipAng_P=new TH1F**    [NumSlicesBGSat];
00110   for(Int_t ck=0;ck<NumSlicesBGSat;ck++)
00111     { 
00112       sprintf(NameHist,"QvsDipAngP%d",ck+1);
00113       QvsDipAngP[ck]=new TProfile(NameHist,"dE/dx",DipAngSlices,DipAngMin,DipAngLimit);
00114       QvsDipAngP[ck]->SetErrorOption("s");
00115       QvsDipAng_P[ck]=new   TH1F* [DipAngSlices];
00116       for(Int_t ida=0;ida<DipAngSlices;ida++)
00117         { 
00118           sprintf(NameHist,"QvsDipAng_P%d_%d",ida+1,ck+1);
00119           QvsDipAng_P[ck][ida]=new TH1F(NameHist,"dE/dx",NumHistBins, 200, 1500);
00120         }
00121     }
00122   dEdxHist=new TH1D("dEdxHist","dEdxHist",NumHistBins, MinHistValue, MaxHistValue);
00123   //stepP=(BetaGammaLogMaxSat-BetaGammaLogMinSat)/NumSlicesBGSat;
00124   stepP=(MomentMaxSat-MomentMinSat)/NumSlicesBGSat;
00125   DipAngStep=(DipAngLimit-DipAngMin)/DipAngSlices;
00126 }

void DedxCalib::checkSelections (   )  [inherited]
 

void DedxCalib::checkSelections (   )  [inherited]
  

00121                                 {
00122   log<<MSG::INFO<<"DedxCalib::checkSelections()...."<<endreq;
00123 }

virtual StatusCode DedxCalib::execute (   )  [virtual, inherited]
 

StatusCode DedxCalib::execute (   )  [virtual, inherited]
  

00075                               {
00076   
00077   this->FillTree();
00078   
00079   this->FillHists();
00080   return StatusCode::SUCCESS;
00081 }

void DedxCalibDipAngle::FillHists (   )  [virtual]
 

Implements DedxCalib.

void DedxCalibDipAngle::FillHists (   )  [virtual]
 

Implements DedxCalib.

00129 {
00130   log<<MSG::INFO<<"DedxCalibDipAngle::FillHists( )..."<<endreq;  
00131   SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00132   if (!eventHeader) {
00133     log << MSG::FATAL << "Could not find Event Header" << endreq;
00134     return;
00135   } 
00136   int event = eventHeader->eventNumber();
00137   int runNO = eventHeader->runNumber();
00138   int typFlag; //data type flag: 0:MC data,   1: 2150,      2:2200
00139   if(runNO<0) typFlag =0;
00140   else if(runNO>=0 && runNO<5459) typFlag =1;
00141   else if(runNO>=5459 && runNO<8093) typFlag =2;
00142   else if(runNO>=8093) typFlag =3; 
00143   
00144   /*for(Int_t ck=0;ck<NumSlicesBGSat;ck++)
00145     {
00146     QvsDipAngP[ck]->Reset();
00147     for(Int_t ida=0;ida<DipAngSlices;ida++)
00148     {
00149     QvsDipAng_P[ck][ida]->Reset();
00150     }
00151     NumEvInHist[ck]=0;
00152     }
00153   */
00154   Identifier mdcid;
00155   SmartDataPtr<RecMdcTrackCol> trkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00156   if (!trkCol) { 
00157     log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
00158     return;
00159   }
00160   int layid, hid, localwid, lr;
00161   double adc_raw, costhe, zhit, driftdist, dd, ph, sinenta, dipan, pathlength, sintheta, rphi_path;
00162   double m_driftd, driftd; 
00163   double m_P, logp;
00164   //NUMTrk = 0;
00165   //NUMHit = 0;
00166   //  static vector<DedxCalibTrk> ex_trks;
00167   log << MSG::INFO << "Begin to check RecMdcTrackCol"<<endreq; 
00168   RecMdcTrackCol::iterator iter_trk = trkCol->begin();
00169   for( ; iter_trk != trkCol->end(); iter_trk++) {
00170     //DedxCalibTrk trk_ex;
00171     //vector<double> phlist;
00172     //vector<int> layerid;
00173     HepVector a = (*iter_trk)->helix();
00174     HepSymMatrix tkErrM = (*iter_trk)->err();
00175     double m_d0E   = tkErrM[0][0];
00176     double m_phi0E = tkErrM[1][1];
00177     double m_cpaE  = tkErrM[2][2];
00178     double m_z0E   = tkErrM[3][3]; 
00179  
00180     HepPoint3D IP(0,0,0);
00181     Dedx_Helix exhel(IP, a);
00182     //QvsDipAngP->Reset();
00183     log << MSG::DEBUG << "hitList of this  track:" << endreq;
00184     HitRefVec gothits = (*iter_trk)->getVecHits();
00185     HitRefVec::iterator it_gothit = gothits.begin();
00186      
00187     //double m_Pt = 1.0/fabs(a[2]);
00188     //m_P = m_Pt*sqrt(1 + (a[4]*a[4]));
00189     float m_Pt = 1.0/fabs( a(3) );
00190     m_P = m_Pt*sqrt(1 + a(5)*a(5));
00191     if(m_P>1.95||m_P<1.8) continue;
00192     logp = m_P;
00193     //logp=log10(m_P/mass_e);
00194     double m_P_rec = (*iter_trk)->p();
00195     //if(logp<BetaGammaLogMaxSat && logp>BetaGammaLogMinSat && NUMTrk<500000){
00196     if(logp<MomentMaxSat && logp>MomentMinSat && NUMTrk<5000000){
00197     
00198     dipan = atan(a[4]);
00199     int Hit =0;
00200     double sindip = sin(dipan);
00201     SINDip[NUMTrk] = sindip;
00202     LOGP[NUMTrk] = logp;
00203     if(fabs(sindip) > DipAngLimit) continue;   
00204     for( ; it_gothit != gothits.end(); it_gothit++) {
00205       mdcid = (*it_gothit)->getMdcId();
00206       layid =  MdcID::layer(mdcid);
00207       localwid = MdcID::wire(mdcid);
00208       adc_raw = (*it_gothit)->getAdc()*1000000;
00209       zhit = (*it_gothit)->getZhit();
00210       lr = (*it_gothit)->getFlagLR();
00211       if(lr == 0 || lr == 2) m_driftd = (*it_gothit)->getDriftDistLeft();
00212       if(lr == 1 ) m_driftd = (*it_gothit)->getDriftDistRight();
00213       driftd = abs(m_driftd);    
00214       dd = fabs((*it_gothit)->getDoca());
00215       sinenta = sin((*it_gothit)->getEntra());  
00216       costhe = cos(M_PI_2-atan(a[4]));
00217       sintheta = sin(M_PI_2-atan(a[4]));      
00218       int ntpFlag = 1;
00219       ph = exsvc->StandardCorrec(typFlag, ntpFlag, runNO, exhel, mdcid, adc_raw, dd, sinenta, zhit, costhe );
00220       pathlength = exsvc->PathL(ntpFlag, exhel, layid, localwid, zhit);
00221       //pathlength = exsvc->PathLCosmic(exhel, layid, localwid, zhit,m_z0E);   
00222       rphi_path = pathlength * sintheta;
00223       if(layid <8) continue;
00224        // {if(adc_raw<MinADCValuecut || adc_raw>MaxADCValuecut || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;}
00225       else if(layid >= 8)
00226         {if(adc_raw<MinADCValuecut || adc_raw>MaxADCValuecut || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut) continue;}      
00227       if(ph >0){
00228          //QRaw[NUMHit] = adc_raw;
00229          //NormQ[NUMHit] = ph;
00230          QRaw[NUMHit] = ph;
00231          NormQ[NUMHit] = adc_raw; 
00232        }
00233       NUMHit++;
00234       Hit++;
00235     }
00236     NUMHits[NUMTrk] = Hit;      
00237     NUMTrk++;  
00238    }
00239   }
00240 }

void DedxCalibDipAngle::FillHistsFromTree (   )  [virtual]
 

Implements DedxCalib.

void DedxCalibDipAngle::FillHistsFromTree (   )  [virtual]
 

Implements DedxCalib.

00595 {
00596   log<<MSG::INFO<<"DedxCalibDipAngle::FillHistsFromTree( )..."<<endreq;  
00597 }

void DedxCalib::FillTestHists (   )  [inherited]
 

void DedxCalib::FillTestHists (   )  [inherited]
  

00117                               {
00118   log<<MSG::INFO<<"DedxCalib::FillTestHists()...."<<endreq;
00119 }

void DedxCalib::FillTree (   )  [inherited]
 

void DedxCalib::FillTree (   )  [inherited]
  

00105                          {
00106   log<<MSG::INFO<<"DedxCalib::FillTree()...."<<endreq;
00107 }

virtual StatusCode DedxCalib::finalize (   )  [virtual, inherited]
 

StatusCode DedxCalib::finalize (   )  [virtual, inherited]
  

00084                                {
00085   
00086   log << MSG::INFO << "DedxCalib finalize() ..." << endreq;
00087   this->AnalyseHists();
00088   
00089   this->WriteCalibdEdxParameters();
00090   this->WriteHists();
00091 
00092   this->FillTestHists();
00093 
00094   
00095   
00096   return StatusCode::SUCCESS;
00097 }

virtual Float_t DedxCalibDipAngle::GetChargeOffCorr (  Int_t  il,
Int_t  iw,
Float_t  ChargeIn,
Float_t  path,
Float_t  doca,
Float_t  EAngle,
Float_t  dipAngle,
Float_t  z
)  [virtual]
 

Implements DedxCalib.

Float_t DedxCalibDipAngle::GetChargeOffCorr (  Int_t  il,
Int_t  iw,
Float_t  ChargeIn,
Float_t  path,
Float_t  doca,
Float_t  EAngle,
Float_t  dipAngle,
Float_t  z
)  [virtual]
 

Implements DedxCalib.

00063 {
00064   log<<MSG::INFO<<"DedxCalibDipAngle::GetChargeOffCorr( )..."<<endreq;
00065   return 1.0;
00066 }

Float_t DedxCalib::GetChargeOnCorr (  Int_t  il,
Int_t  iw,
Float_t  ChargeIn,
Float_t  path,
Float_t  doca,
Float_t  EAngle,
Float_t  dipAngle,
Float_t  z
)  [inherited]
 

Float_t DedxCalib::GetChargeOnCorr (  Int_t  il,
Int_t  iw,
Float_t  ChargeIn,
Float_t  path,
Float_t  doca,
Float_t  EAngle,
Float_t  dipAngle,
Float_t  z
)  [inherited]
  

00127                                                                {
00128   
00129 }

virtual StatusCode DedxCalib::initialize (   )  [virtual, inherited]
 

StatusCode DedxCalib::initialize (   )  [virtual, inherited]
  

00035                                  {
00036   log << MSG::INFO << "DedxCalib initialze() ..." << endreq;
00037   /*StatusCode scint = Service::initialize();
00038   if( scint.isFailure() ) return scint;
00039   IIncidentSvc* incsvc;
00040   scint = service("IncidentSvc", incsvc);
00041   int priority = 100;
00042   if( sc.isSuccess() ){
00043     incsvc -> addListener(this, "BeginEvent", priority);
00044    //incsvc -> addListener(this, "NewRun", priority);
00045   }*/  
00046 
00047   StatusCode sc = service("MdcGeomSvc", geosvc);   
00048   if (sc ==  StatusCode::SUCCESS) {              
00049     log << MSG::INFO <<"MdcGeomSvc is running"<<endl;
00050   } else {
00051     log << MSG::ERROR <<"MdcGeomSvc is failed"<<endl;
00052     return StatusCode::SUCCESS;
00053   }
00054 
00055   StatusCode scex = service("DedxCorrecSvc", exsvc);
00056   if (scex ==  StatusCode::SUCCESS) {              
00057     log << MSG::INFO <<"Hi, DedxCorrectSvc is running"<<endl;
00058   } else {
00059     log << MSG::ERROR <<"DedxCorrectSvc is failed"<<endl;
00060     return StatusCode::SUCCESS;
00061   }
00062 
00063   exsvc->set_flag( calib_flag );
00064   this->checkSelections();
00065   log << MSG::INFO <<"DedxCalib: read correction parameters"<<endreq;
00066   this->ReadCalibdEdxParameters();
00067   this->FillHistsFromTree();
00068   
00069   this->BookHists();
00070   
00071   return StatusCode::SUCCESS;
00072 }

void DedxCalib::ReadBetheBlochParameters (   )  [inherited]
 

void DedxCalib::ReadBetheBlochParameters (   )  [inherited]
 

void DedxCalib::ReadCalibdEdxParameters (   )  [inherited]
 

void DedxCalib::ReadCalibdEdxParameters (   )  [inherited]
  

00109                                         {
00110   log<<MSG::INFO<<"DedxCalib::ReadCalibdEdxParameters()...."<<endreq;
00111 }

void DedxCalibDipAngle::WriteCalibdEdxParameters (   )  [virtual]
 

Reimplemented from DedxCalib.

void DedxCalibDipAngle::WriteCalibdEdxParameters (   )  [virtual]
 

Reimplemented from DedxCalib.

00600 {  
00601   log<<MSG::INFO<<"DedxCalibDipAngle::WriteCalibdEdxParameters( )..."<<endreq;
00602 }

void DedxCalibDipAngle::WriteHists (   )  [virtual]
 

Implements DedxCalib.

void DedxCalibDipAngle::WriteHists (   )  [virtual]
 

Implements DedxCalib.

00069 {
00070   log<<MSG::INFO<<"DedxCalibDipAngle::WriteHists( )..."<<endreq;
00071 
00072   TFile fhist("/ihepbatch/besd09/xcao/calibconst_temp/dug_root/dipcalib.root", "recreate");
00073   TTree* dipar = new TTree("dipcalib", "dipcalib");
00074   if(OldAnalyse==0)
00075   {dipar-> Branch("dipar", &SavePar_new, "dipar[8]/D");}
00076   else {
00077   //cout<<"par old id : "<<SavePar_old[0]<<",         "<<SavePar_old[1]<<",          "<<SavePar_old[2]<<endl;
00078   dipar-> Branch("dipar", &SavePar_old, "dipar[3]/D");}
00079   dipar -> Fill();
00080   dipar -> Write();
00081   
00082   for(Int_t ck=0;ck<NumSlicesBGSat;ck++)
00083     {
00084       QvsDipAngP[ck]->Write();
00085       for(Int_t ida=0;ida<DipAngSlices;ida++)
00086         {
00087           QvsDipAng_P[ck][ida]->Write();
00088         }
00089     }
00090  
00091   delete dipar;
00092   for(Int_t ck=0;ck<NumSlicesBGSat;ck++)
00093     {
00094       delete QvsDipAngP[ck];
00095       for(Int_t ida=0;ida<DipAngSlices;ida++)
00096         {
00097           delete QvsDipAng_P[ck][ida];
00098         }
00099     }
00100 
00101 }

Friends And Related Function Documentation

void fcnQNew (  Int_t &  npar,
Double_t *  gin,
Double_t &  f,
Double_t *  par,
Int_t  iflag
)  [friend]
  

00502                                                                                  {
00503   
00504   Double_t chisq=0.0;
00505   Double_t sinDipAng,p,averdE=0,logp,meandE;
00506   Double_t nentries=0;
00507   Double_t dap[5];
00508   std::vector<Float_t>* vde=0;
00509   Int_t           iP,NumSamples,numh,igate;
00510   Double_t  qtmp,path,normcharge;
00511   Double_t  mult;
00512   for(Int_t ck=0;ck<NumSlicesBGSat;ck++)
00513     {
00514       QvsDipAngP[ck]->Reset();
00515  
00516     }
00517   dap[0]=par[3];
00518   dap[1]=par[4];
00519   dap[2]=par[5];
00520   dap[3]=par[6];
00521   dap[4]=par[7];
00522   long int  nh=0;
00523   for(Int_t i=0;i<NUMTrk;i++)
00524     {
00525       iP=(Int_t)floor((LOGP[i]-MomentMinSat)/stepP);
00526       sinDipAng=fabs(SINDip[i]);
00527       vde=new std::vector<Float_t> (NUMHits[i]);
00528       mult=1./sqrt(sq(sinDipAng)+sq(0.1*par[1]));
00529       
00530       for(Int_t j=0;j<NUMHits[i];j++)
00531         {
00532           (*vde)[j]=QRaw[nh]*(1+par[0]*1e-4*QRaw[nh]*mult+par[2]*1e-09*sq(QRaw[nh]*mult))
00533             *NormQ[nh]/ChebyP4_dip(&sinDipAng,dap);//fix it
00534           nh++;
00535          //cout <<"QRaw["<<nh<<"]="<< QRaw[nh]<<"        NormQ[" <<nh<<"]="<<NormQ[nh]<<" ChebyP4= "<<ChebyP4<<endl;
00536          //cout <<"mult = "<<mult << "         (*vde)[j] = :  " <<(*vde)[j]<< endl; 
00537        }
00538       std::stable_sort(vde->begin(),vde->end());
00539       NumSamples=(Int_t)floor(0.8*NUMHits[i]);
00540       for(Int_t ns=0;ns<NumSamples;ns++) averdE+=(*vde)[ns];
00541       averdE/=NumSamples;
00542       //cout<<"averdE = "<<averdE<<endl;
00543       QvsDipAngP[iP]->Fill(fabs(sinDipAng),averdE);
00544       averdE=0;
00545       delete vde;
00546     }
00547   //log<<MSG::INFO << " nh: " <<nh<< endreq;
00548   //cout << " nh: " <<nh<< endl; 
00549   chisq=0.0;
00550  
00551   Int_t np=0;
00552   for(Int_t k=0;k<NumSlicesBGSat;k++)
00553     {
00554       if(NumEvInHist[k]>40&&QvsDipAngP[k]->GetMean(2)>0&&QvsDipAngP[k]->GetRMS(2)>0)
00555         {
00556           meandE=QvsDipAngP[k]->GetMean(2);
00557           for(Int_t id=0;id<DipAngSlices;id++)
00558             {
00559  
00560               if(QvsDipAngP[k]->GetBinEntries(id)>40)
00561                 {
00562                   chisq+=sq(QvsDipAngP[k]->GetBinContent(id)-meandE)/sq(QvsDipAngP[k]->GetBinError(id));
00563                   np++;
00564                 }
00565             }
00566          // log<<MSG::INFO << " np: "<<np<< endreq ;
00567             cout<< " np: "<<np<< endl;  
00568         }
00569     }
00570   /*log<<MSG::INFO<< "  p0: " <<setprecision(12)<< par[0]
00571                 << "  p1: " << setprecision(12) <<  par[1]
00572                 << "  p2: " << setprecision(14)<< par[2]
00573                 << "  a1: " << setprecision(12)<< par[3]
00574                 << "  a2: " << setprecision(12)<<  par[4]
00575                 << "  a3: " << setprecision(12)<< par[5]
00576                 << "  a4: " << setprecision(12)<< par[6]
00577                 << "  a5: " << setprecision(12)<< par[7]
00578                 << "  chisq: " <<setprecision(12)<< chisq << endreq;
00579   */
00580   
00581            cout << "  p0: " << setprecision(12)<< par[0]
00582                 << "  p1: " << setprecision(12)<< par[1]
00583                 << "  p2: " << setprecision(14)<< par[2]
00584                 << "  a1: " << setprecision(12)<< par[3]
00585                 << "  a2: " << setprecision(12)<< par[4]
00586                 << "  a3: " << setprecision(12)<< par[5]
00587                 << "  a4: " << setprecision(12)<< par[6]
00588                 << "  a5: " << setprecision(12)<< par[7]
00589                 << "  chisq: " <<setprecision(12)<< chisq << endl;
00590 
00591   f=chisq ;
00592 }

void fcnQNew (  Int_t &  npar,
Double_t *  gin,
Double_t &  f,
Double_t *  par,
Int_t  iflag
)  [friend]
  

00502                                                                                  {
00503   
00504   Double_t chisq=0.0;
00505   Double_t sinDipAng,p,averdE=0,logp,meandE;
00506   Double_t nentries=0;
00507   Double_t dap[5];
00508   std::vector<Float_t>* vde=0;
00509   Int_t           iP,NumSamples,numh,igate;
00510   Double_t  qtmp,path,normcharge;
00511   Double_t  mult;
00512   for(Int_t ck=0;ck<NumSlicesBGSat;ck++)
00513     {
00514       QvsDipAngP[ck]->Reset();
00515  
00516     }
00517   dap[0]=par[3];
00518   dap[1]=par[4];
00519   dap[2]=par[5];
00520   dap[3]=par[6];
00521   dap[4]=par[7];
00522   long int  nh=0;
00523   for(Int_t i=0;i<NUMTrk;i++)
00524     {
00525       iP=(Int_t)floor((LOGP[i]-MomentMinSat)/stepP);
00526       sinDipAng=fabs(SINDip[i]);
00527       vde=new std::vector<Float_t> (NUMHits[i]);
00528       mult=1./sqrt(sq(sinDipAng)+sq(0.1*par[1]));
00529       
00530       for(Int_t j=0;j<NUMHits[i];j++)
00531         {
00532           (*vde)[j]=QRaw[nh]*(1+par[0]*1e-4*QRaw[nh]*mult+par[2]*1e-09*sq(QRaw[nh]*mult))
00533             *NormQ[nh]/ChebyP4_dip(&sinDipAng,dap);//fix it
00534           nh++;
00535          //cout <<"QRaw["<<nh<<"]="<< QRaw[nh]<<"        NormQ[" <<nh<<"]="<<NormQ[nh]<<" ChebyP4= "<<ChebyP4<<endl;
00536          //cout <<"mult = "<<mult << "         (*vde)[j] = :  " <<(*vde)[j]<< endl; 
00537        }
00538       std::stable_sort(vde->begin(),vde->end());
00539       NumSamples=(Int_t)floor(0.8*NUMHits[i]);
00540       for(Int_t ns=0;ns<NumSamples;ns++) averdE+=(*vde)[ns];
00541       averdE/=NumSamples;
00542       //cout<<"averdE = "<<averdE<<endl;
00543       QvsDipAngP[iP]->Fill(fabs(sinDipAng),averdE);
00544       averdE=0;
00545       delete vde;
00546     }
00547   //log<<MSG::INFO << " nh: " <<nh<< endreq;
00548   //cout << " nh: " <<nh<< endl; 
00549   chisq=0.0;
00550  
00551   Int_t np=0;
00552   for(Int_t k=0;k<NumSlicesBGSat;k++)
00553     {
00554       if(NumEvInHist[k]>40&&QvsDipAngP[k]->GetMean(2)>0&&QvsDipAngP[k]->GetRMS(2)>0)
00555         {
00556           meandE=QvsDipAngP[k]->GetMean(2);
00557           for(Int_t id=0;id<DipAngSlices;id++)
00558             {
00559  
00560               if(QvsDipAngP[k]->GetBinEntries(id)>40)
00561                 {
00562                   chisq+=sq(QvsDipAngP[k]->GetBinContent(id)-meandE)/sq(QvsDipAngP[k]->GetBinError(id));
00563                   np++;
00564                 }
00565             }
00566          // log<<MSG::INFO << " np: "<<np<< endreq ;
00567             cout<< " np: "<<np<< endl;  
00568         }
00569     }
00570   /*log<<MSG::INFO<< "  p0: " <<setprecision(12)<< par[0]
00571                 << "  p1: " << setprecision(12) <<  par[1]
00572                 << "  p2: " << setprecision(14)<< par[2]
00573                 << "  a1: " << setprecision(12)<< par[3]
00574                 << "  a2: " << setprecision(12)<<  par[4]
00575                 << "  a3: " << setprecision(12)<< par[5]
00576                 << "  a4: " << setprecision(12)<< par[6]
00577                 << "  a5: " << setprecision(12)<< par[7]
00578                 << "  chisq: " <<setprecision(12)<< chisq << endreq;
00579   */
00580   
00581            cout << "  p0: " << setprecision(12)<< par[0]
00582                 << "  p1: " << setprecision(12)<< par[1]
00583                 << "  p2: " << setprecision(14)<< par[2]
00584                 << "  a1: " << setprecision(12)<< par[3]
00585                 << "  a2: " << setprecision(12)<< par[4]
00586                 << "  a3: " << setprecision(12)<< par[5]
00587                 << "  a4: " << setprecision(12)<< par[6]
00588                 << "  a5: " << setprecision(12)<< par[7]
00589                 << "  chisq: " <<setprecision(12)<< chisq << endl;
00590 
00591   f=chisq ;
00592 }

void fcnQOLD (  Int_t &  npar,
Double_t *  gin,
Double_t &  f,
Double_t *  par,
Int_t  iflag
)  [friend]
  

00435                                                                                  {
00436 
00437   Double_t chisq=0.0;
00438   Double_t sinDipAng,p,averdE=0,logp,is,meandE;
00439   Double_t nentries=0;
00440   std::vector<Float_t>* vde=0;
00441   Int_t    iP,NumSamples,numh,igate;
00442   Double_t  qtmp,path;
00443   Double_t  mult;
00444   for(Int_t ck=0;ck<NumSlicesBGSat;ck++)
00445   {
00446       QvsDipAngP[ck]->Reset();
00447   }
00448   int long nh=0;
00449   for(Int_t i=0;i<NUMTrk;i++)
00450     {
00451       iP=(Int_t)floor((LOGP[i]-MomentMinSat)/stepP);
00452       sinDipAng=fabs(SINDip[i]);
00453       vde=new std::vector<Float_t> (NUMHits[i]);
00454       mult=1./sqrt(sq(sinDipAng)+sq(0.1*par[1]));
00455       for(Int_t j=0;j<NUMHits[i];j++)
00456         {
00457           (*vde)[j]=QRaw[nh]*(1+par[0]*1e-4*QRaw[nh]*mult+par[2]*1e-09*sq(QRaw[nh]*mult));
00458           nh++;
00459           //cout <<"mult = "<<mult << "         (*vde)[j] = :  " <<(*vde)[j]<< endl; 
00460        }
00461       std::stable_sort(vde->begin(),vde->end());
00462       NumSamples=(Int_t)floor(0.8*NUMHits[i]);
00463       for(Int_t ns=0;ns<NumSamples;ns++) averdE+=(*vde)[ns];
00464       averdE/=NumSamples;
00465       //cout<<"averdE = "<<averdE<<"      iP= "<<iP<<endl;
00466       QvsDipAngP[iP]->Fill(sinDipAng,averdE,1.);
00467       averdE=0;
00468       delete vde;
00469     }
00470   //log<<MSG::INFO << " nh: " <<nh<< endreq;
00471     chisq=0.0;
00472     for(Int_t k=0;k<NumSlicesBGSat;k++)
00473     {  
00474        //if(NumEvInHist[k]>20&&QvsDipAngP[k]->GetMean(2)>0&&QvsDipAngP[k]->GetRMS(2)>0)
00475        if(QvsDipAngP[k]->GetMean(2)>0&&QvsDipAngP[k]->GetRMS(2)>0)
00476        { 
00477           meandE=QvsDipAngP[k]->GetMean(2);
00478           //log<<MSG::INFO << " Qmean: "<<meandE<<endreq;
00479           for(Int_t id=0;id<DipAngSlices;id++)
00480             { //cout<<"QvsDipAngP[k]->GetBinEntries(id) = "<<QvsDipAngP[k]->GetBinEntries(id)<<endl;
00481               if(QvsDipAngP[k]->GetBinEntries(id)>40)
00482                {
00483                   chisq+=sq(QvsDipAngP[k]->GetBinContent(id)-meandE)/sq(QvsDipAngP[k]->GetBinError(id));
00484  
00485                 }
00486             }
00487         }
00488     }
00489     cout<< " p0: " <<setprecision(12)<< par[0]
00490             << " p1: " <<setprecision(12) <<par[1]
00491             << " p2: " <<setprecision(12) <<par[2]
00492             << " chisq: " <<setprecision(12)<< chisq << endl;
00493 
00494     /*log<<MSG::INFO << " p0: " <<setprecision(12)<< par[0]
00495                  << " p1: " <<setprecision(12) <<par[1]
00496                  << " p2: " <<setprecision(12) <<par[2]
00497                  << " chisq: " <<setprecision(12)<< chisq << endreq;*/
00498     f=chisq ;  
00499 }

void fcnQOLD (  Int_t &  npar,
Double_t *  gin,
Double_t &  f,
Double_t *  par,
Int_t  iflag
)  [friend]
  

00435                                                                                  {
00436 
00437   Double_t chisq=0.0;
00438   Double_t sinDipAng,p,averdE=0,logp,is,meandE;
00439   Double_t nentries=0;
00440   std::vector<Float_t>* vde=0;
00441   Int_t    iP,NumSamples,numh,igate;
00442   Double_t  qtmp,path;
00443   Double_t  mult;
00444   for(Int_t ck=0;ck<NumSlicesBGSat;ck++)
00445   {
00446       QvsDipAngP[ck]->Reset();
00447   }
00448   int long nh=0;
00449   for(Int_t i=0;i<NUMTrk;i++)
00450     {
00451       iP=(Int_t)floor((LOGP[i]-MomentMinSat)/stepP);
00452       sinDipAng=fabs(SINDip[i]);
00453       vde=new std::vector<Float_t> (NUMHits[i]);
00454       mult=1./sqrt(sq(sinDipAng)+sq(0.1*par[1]));
00455       for(Int_t j=0;j<NUMHits[i];j++)
00456         {
00457           (*vde)[j]=QRaw[nh]*(1+par[0]*1e-4*QRaw[nh]*mult+par[2]*1e-09*sq(QRaw[nh]*mult));
00458           nh++;
00459           //cout <<"mult = "<<mult << "         (*vde)[j] = :  " <<(*vde)[j]<< endl; 
00460        }
00461       std::stable_sort(vde->begin(),vde->end());
00462       NumSamples=(Int_t)floor(0.8*NUMHits[i]);
00463       for(Int_t ns=0;ns<NumSamples;ns++) averdE+=(*vde)[ns];
00464       averdE/=NumSamples;
00465       //cout<<"averdE = "<<averdE<<"      iP= "<<iP<<endl;
00466       QvsDipAngP[iP]->Fill(sinDipAng,averdE,1.);
00467       averdE=0;
00468       delete vde;
00469     }
00470   //log<<MSG::INFO << " nh: " <<nh<< endreq;
00471     chisq=0.0;
00472     for(Int_t k=0;k<NumSlicesBGSat;k++)
00473     {  
00474        //if(NumEvInHist[k]>20&&QvsDipAngP[k]->GetMean(2)>0&&QvsDipAngP[k]->GetRMS(2)>0)
00475        if(QvsDipAngP[k]->GetMean(2)>0&&QvsDipAngP[k]->GetRMS(2)>0)
00476        { 
00477           meandE=QvsDipAngP[k]->GetMean(2);
00478           //log<<MSG::INFO << " Qmean: "<<meandE<<endreq;
00479           for(Int_t id=0;id<DipAngSlices;id++)
00480             { //cout<<"QvsDipAngP[k]->GetBinEntries(id) = "<<QvsDipAngP[k]->GetBinEntries(id)<<endl;
00481               if(QvsDipAngP[k]->GetBinEntries(id)>40)
00482                {
00483                   chisq+=sq(QvsDipAngP[k]->GetBinContent(id)-meandE)/sq(QvsDipAngP[k]->GetBinError(id));
00484  
00485                 }
00486             }
00487         }
00488     }
00489     cout<< " p0: " <<setprecision(12)<< par[0]
00490             << " p1: " <<setprecision(12) <<par[1]
00491             << " p2: " <<setprecision(12) <<par[2]
00492             << " chisq: " <<setprecision(12)<< chisq << endl;
00493 
00494     /*log<<MSG::INFO << " p0: " <<setprecision(12)<< par[0]
00495                  << " p1: " <<setprecision(12) <<par[1]
00496                  << " p2: " <<setprecision(12) <<par[2]
00497                  << " chisq: " <<setprecision(12)<< chisq << endreq;*/
00498     f=chisq ;  
00499 }

Member Data Documentation

int DedxCalib::calib_flag [protected, inherited]
 

bool DedxCalib::ddgflag [protected, inherited]
 

IDedxCorrecSvc* DedxCalib::exsvc [protected, inherited]
 

IDedxCorrecSvc* DedxCalib::exsvc [protected, inherited]
 

IMdcGeomSvc* DedxCalib::geosvc [protected, inherited]
 

IMdcGeomSvc* DedxCalib::geosvc [protected, inherited]
 

bool DedxCalib::ggsflag [protected, inherited]
 

bool DedxCalib::layergflag [protected, inherited]
 

MsgStream DedxCalib::log [protected, inherited]
 

TObjArray* DedxCalibDipAngle::m_caliblist [private]
 

TObjArray* DedxCalibDipAngle::m_caliblist [private]
 

std::string DedxCalib::m_constrootfile [protected, inherited]
 

int DedxCalib::m_correc_flag [protected, inherited]
 

std::string DedxCalib::m_inputfile [protected, inherited]
 

double DedxCalibDipAngle::m_par
 

int DedxCalib::m_par_flag [protected, inherited]
 

int DedxCalib::m_phshape_flag [protected, inherited]
 

TRandom* DedxCalibDipAngle::m_rand [private]
 

TRandom* DedxCalibDipAngle::m_rand [private]
 

std::string DedxCalib::m_rootfile [protected, inherited]
 

TObjArray* DedxCalibDipAngle::m_satlist [private]
 

TObjArray* DedxCalibDipAngle::m_satlist [private]
 

bool DedxCalib::wiregflag [protected, inherited]
 

bool DedxCalib::zdepflag [protected, inherited]
 

The documentation for this class was generated from the following files:
Generated on Wed Feb 2 15:57:05 2011 for BOSS6.5.5 by  doxygen 1.3.9.1