BOSS6.5.5: DedxCalibLayerGain Class Reference

00028   : DedxCalib(name, pSvcLocator)
00029 {
00030   log<<MSG::INFO<<"DedxCalibLayerGain::DedxCalibLayerGain( )..."<<endreq;
00031 }

00349 { 
00350   log<<MSG::INFO<<"DedxCalibLayerGain::AnalyseHists( )..."<<endreq;
00351   TF1* func;
00352   if(m_phshape_flag==0) {
00353     func = new TF1("mylan",mylan,1000,4000,4);
00354     func->SetParameters(10000.0,2010.1,300,-1.0);
00355   } else if(m_phshape_flag==1) { 
00356     func = new TF1("landaun",landaun,1000,4000,3);
00357     func->SetParameters(2000.0,300.0);
00358   } else if(m_phshape_flag==2) { 
00359     func = new TF1("land",Landau,1000,4000,2);
00360     func->SetParameters(2010.1,300);
00361   } else if(m_phshape_flag==3) { 
00362     func = new TF1("vavi",Vavilov,1000,4000,2);
00363     func->SetParameters(3.1,0.4);
00364   } else if(m_phshape_flag==4) { 
00365     func = new TF1("asym",AsymGauss,1000,4000,4);
00366     func->SetParameters(10000.0,2000.0,300.0,100.);
00367   } 
00368 
00369   func->SetLineWidth(2.1);
00370   func->SetLineColor(2);   
00371   Double_t entry, mean, rms;
00372   Double_t binmax;
00373   Double_t lp[3];
00374   for(int lay=0; lay<43; lay++){
00375     //cout<<"m_phshape_flag = "<<m_phshape_flag<<endl;
00376     entry = m_hist[lay]->GetEntries();
00377     mean = m_hist[lay]->GetMean();
00378     rms = m_hist[lay]->GetRMS();
00379     binmax = m_hist[lay]->GetMaximumBin();
00380     lp[1] = m_hist[lay]->GetBinCenter(binmax);
00381     lp[2] = 200;
00382     lp[0] = (m_hist[lay]->GetBinContent(binmax)+m_hist[lay]->GetBinContent(binmax-1)+m_hist[lay]->GetBinContent(binmax+1))/3;
00383                    
00384     if(m_phshape_flag==0) {
00385       func->SetParameters(entry, mean, rms, -0.6);
00386     } else if(m_phshape_flag==1) {     
00387       //func->SetParameters(entry, 0.7*mean, rms );
00388       func->SetParameters(lp[0], lp[1], lp[2] );
00389     } else if(m_phshape_flag==2) {     
00390       func->SetParameters(0.7*mean, rms );
00391     } else if(m_phshape_flag==3) { 
00392       func->SetParameters(0.05, 0.6);
00393     } else if(m_phshape_flag==4) {     
00394       func->SetParameters(entry, mean, rms, 1.0 );
00395     }
00396     m_hist[lay]->Fit(func, "Q" );    
00397   }                           
00398 }

00225 {
00226   log<<MSG::INFO<<"DedxCalibLayerGain::BookHists( )..."<<endreq;  
00227   log<<MSG::DEBUG<<" bookhist...1" <<endreq;                   
00228   m_layerglist = new TObjArray(0);
00229   log<<MSG::DEBUG<<" bookhist...2" <<endreq;                   
00230   m_rand = new TRandom();
00231   std::string hlabel;            
00232   log<<MSG::DEBUG<<" bookhist...3" <<endreq;
00233 
00234   for(int jj=0; jj<43; jj++) {
00235     std::ostringstream strout;                                                 
00236     strout << "Layer_Gain_" << jj;
00237     hlabel = strout.str();
00238     if(jj<8){
00239        TH1F *layerg = new TH1F(hlabel.c_str(), hlabel.c_str(), NumHistBins, Iner_MinHistValue, Iner_MaxHistValue);
00240        layerg->GetXaxis()->SetTitle("dE/dx(eV/1.5cm)");
00241        layerg->GetYaxis()->SetTitle("entries/(40eV/1.5cm)");
00242        m_hist.push_back(layerg);
00243        m_layerglist -> Add(layerg);
00244     }
00245     else if(jj>=8){
00246        TH1F *layerg = new TH1F(hlabel.c_str(), hlabel.c_str(), NumHistBins, MinHistValue, MaxHistValue);  
00247        //    layerg->Print("base");
00248        layerg->GetXaxis()->SetTitle("dE/dx(eV/1.5cm)");
00249        layerg->GetYaxis()->SetTitle("entries/(40eV/1.5cm)");
00250        m_hist.push_back(layerg);                                              
00251        m_layerglist -> Add(layerg);
00252     }                             
00253   }
00254   log<<MSG::DEBUG<<" m_hist size..." <<m_hist.size()<<endreq;  
00255 }

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

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

00258 {
00259   log<<MSG::INFO<<"DedxCalibLayerGain::FillHists( )..."<<endreq;  
00260   SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00261   if (!eventHeader) {
00262     log << MSG::FATAL << "Could not find Event Header" << endreq;
00263     return;
00264   } 
00265   int event = eventHeader->eventNumber();
00266   int runNO = eventHeader->runNumber();
00267   int typFlag; //data type flag: 0:MC data,   1: 2150,      2:2200
00268   if(runNO<0) typFlag =0;
00269   else if(runNO>=0 && runNO<5459) typFlag =1;
00270   else if(runNO>=5459 && runNO<8093) typFlag =2;
00271   else if(runNO>=8093 && runNO<9815) typFlag =3;
00272   else if(runNO>=9815 ) typFlag =4;
00273  
00274   double tes;
00275   int esTimeflag = -1;
00276   SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00277   if( ! estimeCol) {
00278           log << MSG::WARNING << "Could not find EvTimeCol" << endreq;
00279   }else{
00280           RecEsTimeCol::iterator iter_evt = estimeCol->begin();
00281           for(; iter_evt!=estimeCol->end(); iter_evt++){
00282                   tes =  (*iter_evt)->getTest();
00283                   esTimeflag = (*iter_evt)->getStat();
00284           }
00285   }
00286   if(typFlag ==2) {if( tes>1000 ) return;}
00287   else if(typFlag ==3 ){if (tes>700 ) return;}
00288   else if(typFlag ==4 ){if (tes>1400 ) return;}
00289 
00290   Identifier mdcid;
00291   SmartDataPtr<RecMdcTrackCol> trkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00292   if (!trkCol) { 
00293      log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
00294      return;
00295   }
00296   int lr, layid, hid,localwid,w0id,wid;
00297   double adc_raw, costhe, zhit, ddl, ddr, dd, ph, sinenta, pathlength, sintheta, rphi_path;
00298   double m_driftd, driftd;
00299   log << MSG::INFO << "Begin to check RecMdcTrackCol"<<endreq; 
00300   RecMdcTrackCol::iterator iter_trk = trkCol->begin();
00301    for( ; iter_trk != trkCol->end(); iter_trk++) {
00302       HepVector a = (*iter_trk)->helix();
00303       HepSymMatrix tkErrM = (*iter_trk)->err();
00304       double m_d0E   = tkErrM[0][0];
00305       double m_phi0E = tkErrM[1][1];
00306       double m_cpaE  = tkErrM[2][2];
00307       double m_z0E   = tkErrM[3][3];  
00308       HepPoint3D IP(0,0,0);
00309       Dedx_Helix exhel(IP, a);
00310       log << MSG::DEBUG << "hitList of this  track:" << endreq;
00311       HitRefVec gothits = (*iter_trk)->getVecHits();
00312       HitRefVec::iterator it_gothit = gothits.begin();
00313       for( ; it_gothit != gothits.end(); it_gothit++) {
00314         mdcid = (*it_gothit)->getMdcId();
00315         layid =  MdcID::layer(mdcid);
00316         localwid = MdcID::wire(mdcid);
00317         w0id = geosvc->Layer(layid)->Wirst();
00318         wid = w0id + localwid;  
00319         adc_raw = (*it_gothit)->getAdc();
00320         zhit = (*it_gothit)->getZhit();
00321         lr = (*it_gothit)->getFlagLR();
00322         if(lr == 0 || lr == 2) m_driftd = (*it_gothit)->getDriftDistLeft();
00323         if(lr == 1 ) m_driftd = (*it_gothit)->getDriftDistRight();
00324         driftd = abs(m_driftd);
00325         dd = fabs((*it_gothit)->getDoca());
00326         sinenta = sin((*it_gothit)->getEntra());        
00327         costhe = cos(M_PI_2-atan(a[4]));
00328         sintheta = sin(M_PI_2-atan(a[4]));
00329         int ntpFlag = 1;
00330         ph = exsvc->StandardCorrec(typFlag, ntpFlag, runNO, exhel, mdcid, adc_raw, dd, sinenta, zhit, costhe );
00331         pathlength = exsvc->PathL(ntpFlag, exhel, layid, localwid, zhit);
00332         //ph = exsvc->StandardCorrec( exhel, mdcid, adc_raw, dd, sinenta, zhit, costhe ); 
00333         //pathlength = exsvc->PathL(exhel, layid, localwid, zhit);
00334         //pathlength = exsvc->PathLCosmic(exhel, layid, localwid, zhit,m_z0E);
00335         rphi_path = pathlength * sintheta;
00336         //cout<<layid<<"     "<<localwid<<"     "<<ph<<endl;
00337         if(layid <8)
00338         {if( adc_raw<MinADCValuecut || adc_raw>MaxADCValuecut || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;}
00339         else if(layid >= 8)
00340         {if(adc_raw<MinADCValuecut || adc_raw>MaxADCValuecut || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut) continue;}       
00341                 
00342         m_hist[layid]->Fill(ph);                                 
00343       } 
00344    }
00345 }

00401 {
00402   log<<MSG::INFO<<"DedxCalibLayerGain::FillHistsFromTree( )..."<<endreq;  
00403 }

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

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

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 }

00036 {
00037   log<<MSG::INFO<<"DedxCalibLayerGain::GetChargeOffCorr( )..."<<endreq;
00038   return 1.0;
00039 }

00127                                                                {
00128   
00129 }

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 }

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

00406 {  
00407   log<<MSG::INFO<<"DedxCalibLayerGain::WriteCalibdEdxParameters( )..."<<endreq;
00408 }

00043 {
00044   log<<MSG::INFO<<"DedxCalibLayerGain::WriteHists( )..."<<endreq;
00045   //double mean[43],mean0,layNO[43],meanlay[43];
00046   double layermean[43],layermean0,layergain[43],layerchisq[43],meanlayer[43],rmslayer[43],sigmalayer[43];
00047   double layerentries[43];
00048   double layNO[43];
00049       
00050   if(m_phshape_flag==0) {
00051     layermean0 = m_hist[37] -> GetFunction("mylan") -> GetParameter(1); 
00052   } else   if(m_phshape_flag==1) {
00053     layermean0 = m_hist[37] -> GetFunction("landaun") -> GetParameter(1); 
00054   } else   if(m_phshape_flag==2) {
00055     layermean0 = m_hist[37] -> GetFunction("land") -> GetParameter(1);
00056   } else   if(m_phshape_flag==3) {
00057     layermean0 = m_hist[37] -> GetFunction("vavi") -> GetParameter(1); 
00058   } else   if(m_phshape_flag==4) {
00059     layermean0 = m_hist[37] -> GetFunction("asym") -> GetParameter(1); 
00060   }
00061 
00062   for(int lay=0; lay<43; lay++) {
00063      layerentries[lay] = m_hist[lay]->GetEntries();
00064      meanlayer[lay] = m_hist[lay]->GetMean();
00065      rmslayer[lay] = m_hist[lay]->GetRMS();
00066      layNO[lay] = lay;
00067     if(m_phshape_flag==0) {
00068       layermean[lay]  = m_hist[lay] -> GetFunction("mylan") -> GetParameter(1); 
00069       layerchisq[lay] = m_hist[lay] -> GetFunction("mylan") -> GetParameter(2);
00070       sigmalayer[lay] = m_hist[lay] -> GetFunction("mylan") -> GetParError(1);
00071     } else if(m_phshape_flag==1) {
00072       layermean[lay]  = m_hist[lay] -> GetFunction("landaun") -> GetParameter(1); 
00073       layerchisq[lay] = m_hist[lay] -> GetFunction("landaun") -> GetParameter(2);
00074       sigmalayer[lay] = m_hist[lay] -> GetFunction("landaun") -> GetParError(1);
00075     } else if(m_phshape_flag==2) {
00076       layermean[lay]  = m_hist[lay] -> GetFunction("land") -> GetParameter(1); 
00077     } else if(m_phshape_flag==3) {
00078       layermean[lay]  = m_hist[lay] -> GetFunction("vavi") -> GetParameter(1); 
00079     } else if(m_phshape_flag==4) {
00080       layermean[lay]  = m_hist[lay] -> GetFunction("asym") -> GetParameter(1); 
00081     }
00082     //cout<<"lay = "<<layNO[lay] <<"       layermean = "<<layermean[lay] <<endl;
00083     /*if(layermean0) {
00084       layergain[lay] = layermean[lay]/layermean0;
00085     } else {
00086       layergain[lay] = layermean[lay]/1000.0;
00087     }*/
00088     //layergain[lay] = layermean[lay]/NormalMean;
00089     log<<MSG::INFO<<"mean["<<lay<<"] = "<<m_hist[lay]->GetMean()
00090        <<"fitted mean:"<<layermean[lay]<<endreq;           
00091   }
00092   
00093   double total_mean = 0;
00094   int lay_normal_N =0;
00095   for(int j=8; j<43; j++){
00096       //int wireNO_temp = Normal_wire[j];
00097       if(layermean[j]>0&&sigmalayer[j]<10&&sigmalayer[j]>0&&layerchisq[j]>0&&layerchisq[j]<200&&layerentries[j]>1500)
00098       {
00099        total_mean += layermean[j];
00100        lay_normal_N++;
00101       }
00102   }
00103   double Normal_mean = total_mean/lay_normal_N;
00104   cout<<"Normal_mean : "<<Normal_mean<<endl;
00105   for(int lay=0; lay<43; lay++) {
00106    layergain[lay] = layermean[lay]/Normal_mean;
00107    cout<<"layergain["<<lay<<"] : "<<layergain[lay]<<endl;
00108   }
00109     
00110   TCanvas *ex= new TCanvas("excanvas","dE/dx vs layer",200,10,600,400);
00111   ex->SetFillColor(42);
00112   ex->SetGrid(); 
00113   ex->SetLogx(0);
00114   ex->SetLogy(0);
00115 
00116   TGraph * tc = new TGraph(43,layNO,layermean);
00117   tc->SetMarkerStyle(29);
00118   tc->SetTitle("dE/dx vs layer");
00119   //tc->SetMaximum(800.0);
00120   //tc->SetMinimum(400.0);
00121   tc->SetLineColor(2);
00122   tc->SetLineWidth(4);
00123   tc->SetMarkerColor(4);
00124   tc->GetXaxis()->SetTitle("layer");
00125   tc->GetYaxis()->SetTitle("dE/dx");
00126   tc->GetYaxis()->CenterTitle(1);
00127   tc->GetXaxis()->CenterTitle(1);
00128   tc->Draw("PA");
00129 
00130   ex->Update();
00131   //  ex->GetFrame()->SetFillColor(21);
00132   //  ex->GetFrame()->SetBorderSize(12);
00133   ex->Modified();
00134   
00135   if(m_rootfile=="no rootfile")
00136   {
00137      log << MSG::ERROR<<"no sepcified calibration file path of wiregcalib "<<endreq;
00138      //return StatusCode::FAILURE;
00139      return;
00140   }
00141   std::string rootpath_layer = m_rootfile;
00142   std::string hlabel_rootpath_layer;
00143   std::ostringstream strout;
00144   strout << rootpath_layer <<"/layergcalib.root";
00145   hlabel_rootpath_layer = strout.str();
00146   std::cout<<"string layer :    "<<hlabel_rootpath_layer.c_str()<<endl;
00147  
00148   TFile fhist(hlabel_rootpath_layer.c_str(), "recreate");
00149 
00150   TTree* layer_tree = new TTree("layergcalib", "layergcalib");
00151   layer_tree -> Branch("layerg", &layergain, "layerg[43]/D");
00152   layer_tree -> Branch("layerentries", &layerentries, "layerentries[43]/D");
00153   layer_tree -> Branch("meanlayer", &meanlayer, "meanlayer[43]/D");
00154   layer_tree -> Branch("rmslayer", &rmslayer, "rmslayer[43]/D");
00155   layer_tree -> Branch("layerchisq", &layerchisq, "layerchisq[43]/D");
00156   layer_tree -> Branch("layermean", &layermean, "layermean[43]/D");
00157 
00158   //TFile fhist("/ihepbatch/d09/xcao/calibconst_634/layergcalib.root", "recreate");
00159   layer_tree -> Fill();  
00160   layer_tree -> Write();
00161   m_layerglist -> Write();    
00162   ex -> Write();
00163 
00164   fhist.Close();        
00165   for(int i=0; i<43; i++) delete m_hist[i];     
00166   delete m_layerglist;                                                           
00167   //delete layer_tree;
00168   delete m_rand;             
00169 }


Public Member Functions
void	AnalyseHists ()
void	AnalyseHists ()
void	BookHists ()
void	BookHists ()
void	checkSelections ()
void	checkSelections ()
	DedxCalibLayerGain (const std::string &name, ISvcLocator *pSvcLocator)
	DedxCalibLayerGain (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 ()
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
vector< TH1F * >	m_hist
vector< TH1F * >	m_hist
TObjArray *	m_layerglist
TObjArray *	m_layerglist
TRandom *	m_rand
TRandom *	m_rand

DedxCalibLayerGain Class Reference

Public Member Functions

Protected Attributes

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation