BOSS_7.0.2: /home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Mdc/DedxCalibAlg/DedxCalibAlg-00-01-15/src/DedxCalibMomentum.cxx Source File

00001 #include "GaudiKernel/MsgStream.h"
00002 
00003 #include <sstream>
00004 #include <string>
00005 #include "cmath"
00006 #include "TTree.h" 
00007 #include "TH1F.h"
00008 #include "TF1.h" 
00009 #include "TFile.h"
00010 #include "TStyle.h"
00011 #include "TCanvas.h"
00012 
00013 #include "DedxCalibAlg/DedxCalibMomentum.h"
00014 
00015 #define Size  700000
00016 
00017 using namespace std;
00018 
00019 const double pMin=0.2;
00020 const double pMax=2.1;
00021 const int nbins = 50;
00022 double bin_step = (pMax-pMin)/nbins;
00023 const double chihist_min=-10;
00024 const double chihist_max=10;
00025 const double dedxhist_min=350;
00026 const double dedxhist_max=750;
00027 const double hits_min=0;
00028 const double hits_max=35;
00029 
00030 DedxCalibMomentum::DedxCalibMomentum(const std::string& name, ISvcLocator* pSvcLocator): DedxCalib(name, pSvcLocator){}
00031 
00032 void DedxCalibMomentum::BookHists()
00033 {
00034     MsgStream log(msgSvc(), name());
00035     log<<MSG::INFO<<"DedxCalibMomentum::BookHists()"<<endreq;
00036 
00037     ReadRecFileList(); 
00038 
00039     m_chi = new TH1F*[nbins];
00040     m_pos_chi = new TH1F*[nbins];
00041     m_neg_chi = new TH1F*[nbins];
00042     m_dedx = new TH1F*[nbins];
00043     m_pos_dedx = new TH1F*[nbins];
00044     m_neg_dedx = new TH1F*[nbins];
00045     m_hits = new TH1F*[nbins];
00046 
00047     stringstream hlabel;
00048     for(int i=0;i<nbins;i++)
00049     {
00050         hlabel.str("");
00051         hlabel<<"chi"<<i;
00052         m_chi[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, chihist_min, chihist_max); 
00053         hlabel.str("");
00054         hlabel<<"pos_chi"<<i;
00055         m_pos_chi[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, chihist_min, chihist_max);
00056         hlabel.str("");
00057         hlabel<<"neg_chi"<<i;
00058         m_neg_chi[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, chihist_min, chihist_max);
00059         hlabel.str("");
00060         hlabel<<"dedx"<<i;
00061         m_dedx[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, dedxhist_min, dedxhist_max);
00062         hlabel.str(""); 
00063         hlabel<<"pos_dedx"<<i; 
00064         m_pos_dedx[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, dedxhist_min, dedxhist_max);
00065         hlabel.str("");
00066         hlabel<<"neg_dedx"<<i;  
00067         m_neg_dedx[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, dedxhist_min, dedxhist_max);
00068         hlabel.str("");               
00069         hlabel<<"hits"<<i; 
00070         m_hits[i] =  new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, hits_min, hits_max); 
00071     }
00072 
00073     Vec_dedx.clear();
00074     Vec_ptrk.clear();
00075 }
00076 
00077 void DedxCalibMomentum::FillHists()
00078 {
00079     MsgStream log(msgSvc(), name());
00080     log<<MSG::INFO<<"DedxCalibMomentum::FillHists()"<<endreq;
00081 
00082     TFile* f;
00083     TTree* n103;
00084     string runlist;
00085 
00086     int ndedxhit=0;
00087     double dEdx[100]={0},pathlength[100]={0},wid[100]={0},layid[100]={0},dd_in[100]={0},eangle[100]={0},zhit[100]={0};
00088     double dedx=0;
00089     float runNO=0,runFlag=0,costheta=0,tes=0,charge=0,recalg=0,ptrk=0,chidedx=0;
00090     int usedhit=0;
00091     vector<double> phlist;
00092     cout<<"m_recFileVector.size()= "<<m_recFileVector.size()<<endl;
00093     for(int i=0; i<m_recFileVector.size(); i++)
00094     {
00095         runlist = m_recFileVector[i];
00096         f = new TFile(runlist.c_str());
00097         n103 = (TTree*)f->Get("n103");
00098         n103-> SetBranchAddress("ndedxhit", &ndedxhit);
00099         n103-> SetBranchAddress("dEdx_hit",dEdx);
00100         n103-> SetBranchAddress("pathlength_hit",pathlength);
00101         n103-> SetBranchAddress("wid_hit",wid);
00102         n103-> SetBranchAddress("layid_hit",layid);
00103         n103-> SetBranchAddress("dd_in_hit",dd_in);
00104         n103-> SetBranchAddress("eangle_hit",eangle);
00105         n103-> SetBranchAddress("zhit_hit",zhit);
00106         n103-> SetBranchAddress("runNO",&runNO);
00107         n103-> SetBranchAddress("runFlag",&runFlag);
00108         n103-> SetBranchAddress("costheta",&costheta);
00109         n103-> SetBranchAddress("t0",&tes);
00110         n103-> SetBranchAddress("charge",&charge);
00111         n103-> SetBranchAddress("recalg",&recalg);
00112         n103-> SetBranchAddress("ndedxhit",&ndedxhit);
00113         n103-> SetBranchAddress("ptrk",&ptrk);
00114         n103-> SetBranchAddress("chidedx_e",&chidedx);
00115         for(int j=0;j<n103->GetEntries();j++)
00116         {
00117             phlist.clear();
00118             n103->GetEntry(j);
00119             if(ptrk>pMax || ptrk<pMin) continue;
00120             if(tes>1400) continue;
00121             for(int k=0;k<ndedxhit;k++)
00122             {
00123                 dEdx[k] = exsvc->StandardHitCorrec(0,runFlag,2,runNO,pathlength[k],wid[k],layid[k],dEdx[k],dd_in[k],eangle[k],costheta);
00124                 phlist.push_back(dEdx[k]);
00125             }
00126             dedx = cal_dedx_bitrunc(truncate, phlist);
00127             dedx = exsvc->StandardTrackCorrec(0, runFlag, 2, runNO, dedx, costheta, tes);
00128             int ipBin = (Int_t)floor((ptrk-pMin)/bin_step);
00129             m_dedx[ipBin]->Fill(dedx);
00130             if(charge>0) m_pos_dedx[ipBin]->Fill(dedx);
00131             if(charge<0) m_neg_dedx[ipBin]->Fill(dedx);
00132 
00133             usedhit = ndedxhit*truncate;
00134             set_dEdx(1,dedx,recalg,runFlag,usedhit,ptrk,acos(costheta),1.5,vFlag,tes);
00135             double chi = m_chi_ex[ParticleFlag];
00136             //cout<<"chidedx= "<<chidedx<<" particleType= "<<ParticleFlag<<"  new chi= "<<chi<<endl;
00137             m_chi[ipBin]->Fill(chi);
00138             if(charge>0) m_pos_chi[ipBin]->Fill(chi);
00139             if(charge<0) m_neg_chi[ipBin]->Fill(chi);
00140             m_hits[ipBin]->Fill(usedhit);
00141 
00142             Vec_dedx.push_back(dedx);
00143             Vec_ptrk.push_back(ptrk);
00144         }
00145     }
00146 }
00147 
00148 void DedxCalibMomentum::AnalyseHists()
00149 {
00150     MsgStream log(msgSvc(), name());
00151     log<<MSG::INFO<<"DedxCalibMomentum::AnalyseHists()"<<endreq;
00152 
00153     gStyle->SetOptFit(1111);
00154     for(int i=0; i<nbins; i++)
00155     {
00156         m_dedx[i]->Fit("gaus", "Q" );
00157         m_pos_dedx[i]->Fit("gaus", "Q" );
00158         m_neg_dedx[i]->Fit("gaus", "Q" );
00159         m_chi[i]->Fit("gaus", "Q" );
00160         m_pos_chi[i]->Fit("gaus", "Q" );
00161         m_neg_chi[i]->Fit("gaus", "Q" );
00162     }
00163 }
00164 
00165 void DedxCalibMomentum::WriteHists()
00166 {
00167     MsgStream log(msgSvc(), name());
00168     log<<MSG::INFO<<"DedxCalibMomentum::WriteHists()"<< endreq;
00169 
00170     double chientryNo[nbins]={0},chimean[nbins]={0},chimeanerr[nbins]={0},chisigma[nbins]={0},chisq[nbins]={0};
00171     double pos_chientryNo[nbins]={0},pos_chimean[nbins]={0},pos_chimeanerr[nbins]={0},pos_chisigma[nbins]={0},pos_chisq[nbins]={0};
00172     double neg_chientryNo[nbins]={0},neg_chimean[nbins]={0},neg_chimeanerr[nbins]={0},neg_chisigma[nbins]={0},neg_chisq[nbins]={0};
00173     double fitentryNo[nbins]={0},fitmean[nbins]={0},fitmeanerr[nbins]={0},fitsigma[nbins]={0},fitchisq[nbins]={0};
00174     double pos_fitentryNo[nbins]={0},pos_fitmean[nbins]={0},pos_fitmeanerr[nbins]={0},pos_fitsigma[nbins]={0},pos_fitchisq[nbins]={0};
00175     double neg_fitentryNo[nbins]={0},neg_fitmean[nbins]={0},neg_fitmeanerr[nbins]={0},neg_fitsigma[nbins]={0},neg_fitchisq[nbins]={0};
00176     double hits_mean[nbins]={0},hits_sigma[nbins]={0};
00177     double pBin[nbins]={0};
00178 
00179     for(int i=0;i<nbins;i++)
00180     {
00181         pBin[i] = (i+0.5)*bin_step + pMin;
00182 
00183         chientryNo[i] = m_chi[i]->GetEntries();
00184         if(chientryNo[i]>100)
00185         {
00186             chimean[i] = m_chi[i]->GetFunction("gaus")->GetParameter(1);
00187             chimeanerr[i] = m_chi[i]->GetFunction("gaus")->GetParError(1);
00188             chisigma[i] = m_chi[i]->GetFunction("gaus")->GetParameter(2);
00189             chisq[i] = (m_chi[i]->GetFunction("gaus")->GetChisquare())/(m_chi[i]->GetFunction("gaus")->GetNDF());
00190         }
00191         pos_chientryNo[i] = m_pos_chi[i]->GetEntries();
00192         if(pos_chientryNo[i]>100)
00193         {    
00194             pos_chimean[i] = m_pos_chi[i]->GetFunction("gaus")->GetParameter(1);
00195             pos_chimeanerr[i] = m_pos_chi[i]->GetFunction("gaus")->GetParError(1);
00196             pos_chisigma[i] = m_pos_chi[i]->GetFunction("gaus")->GetParameter(2);
00197             pos_chisq[i] = (m_pos_chi[i]->GetFunction("gaus")->GetChisquare())/(m_pos_chi[i]->GetFunction("gaus")->GetNDF());
00198         }
00199         neg_chientryNo[i] = m_neg_chi[i]->GetEntries();
00200         if(neg_chientryNo[i]>100)
00201         {
00202             neg_chimean[i] = m_neg_chi[i]->GetFunction("gaus")->GetParameter(1);
00203             neg_chimeanerr[i] = m_neg_chi[i]->GetFunction("gaus")->GetParError(1);
00204             neg_chisigma[i] = m_neg_chi[i]->GetFunction("gaus")->GetParameter(2);
00205             neg_chisq[i] = (m_neg_chi[i]->GetFunction("gaus")->GetChisquare())/(m_neg_chi[i]->GetFunction("gaus")->GetNDF()); 
00206         }
00207 
00208         fitentryNo[i] = m_dedx[i]->GetEntries();
00209         if(fitentryNo[i]>100)
00210         {
00211             fitmean[i] = m_dedx[i]->GetFunction("gaus")->GetParameter(1);
00212             fitmeanerr[i] = m_dedx[i]->GetFunction("gaus")->GetParError(1);
00213             fitsigma[i] = m_dedx[i]->GetFunction("gaus")->GetParameter(2);
00214             fitchisq[i] = (m_dedx[i]->GetFunction("gaus")->GetChisquare())/(m_dedx[i]->GetFunction("gaus")->GetNDF());
00215         }
00216         pos_fitentryNo[i] = m_pos_dedx[i]->GetEntries(); 
00217         if(pos_fitentryNo[i]>100)
00218         {
00219             pos_fitmean[i] = m_pos_dedx[i]->GetFunction("gaus")->GetParameter(1);
00220             pos_fitmeanerr[i] = m_pos_dedx[i]->GetFunction("gaus")->GetParError(1);
00221             pos_fitsigma[i] = m_pos_dedx[i]->GetFunction("gaus")->GetParameter(2);     
00222             pos_fitchisq[i] = (m_pos_dedx[i]->GetFunction("gaus")->GetChisquare())/(m_pos_dedx[i]->GetFunction("gaus")->GetNDF());
00223         }   
00224         neg_fitentryNo[i] = m_neg_dedx[i]->GetEntries();
00225         if(neg_fitentryNo[i]>100)
00226         {
00227             neg_fitmean[i] = m_neg_dedx[i]->GetFunction("gaus")->GetParameter(1);
00228             neg_fitmeanerr[i] = m_neg_dedx[i]->GetFunction("gaus")->GetParError(1);
00229             neg_fitsigma[i] = m_neg_dedx[i]->GetFunction("gaus")->GetParameter(2);
00230             neg_fitchisq[i] = (m_neg_dedx[i]->GetFunction("gaus")->GetChisquare())/(m_neg_dedx[i]->GetFunction("gaus")->GetNDF());
00231         } 
00232 
00233         hits_mean[i] = m_hits[i]->GetMean();
00234         hits_sigma[i] = m_hits[i]->GetRMS();
00235     }
00236 
00237     double dedx1[Size] = {0};
00238     double ptrk1[Size] = {0};
00239     cout << "Vec_dedx.size() = " << Vec_dedx.size() << endl;
00240     for(int i=0;i<Size && i< Vec_dedx.size();i++)
00241     {       
00242         dedx1[i] = Vec_dedx[i];
00243         ptrk1[i] = Vec_ptrk[i];
00244         //cout<<"i= "<<i<<"dedx= "<<dedx1[i]<<" ptrk= "<<ptrk1[i]<<endl;
00245     } 
00246 
00247     log<<MSG::INFO<<"begin generating root file!!! "<<endreq;    
00248     TFile* f = new TFile(m_rootfile.c_str(), "RECREATE");
00249     for(int i=0;i<nbins;i++)
00250     {
00251         m_chi[i]->Write();
00252         m_pos_chi[i]->Write();
00253         m_neg_chi[i]->Write();
00254         m_dedx[i]->Write();
00255         m_pos_dedx[i]->Write();
00256         m_neg_dedx[i]->Write();
00257         m_hits[i]->Write();
00258     }
00259 
00260     TTree* momInfor = new TTree("momInfor","momInfor");
00261     momInfor->Branch("chientryNo",chientryNo,"chientryNo[50]/D");
00262     momInfor->Branch("chimean",chimean,"chimean[50]/D");
00263     momInfor->Branch("chimeanerr",chimeanerr,"chimeanerr[50]/D");
00264     momInfor->Branch("chisigma",chisigma,"chisigma[50]/D");
00265     momInfor->Branch("chisq",chisq,"chisq[50]/D");
00266     momInfor->Branch("pos_chientryNo",pos_chientryNo,"pos_chientryNo[50]/D");
00267     momInfor->Branch("pos_chimean",pos_chimean,"pos_chimean[50]/D");
00268     momInfor->Branch("pos_chimeanerr",pos_chimeanerr,"pos_chimeanerr[50]/D");
00269     momInfor->Branch("pos_chisigma",pos_chisigma,"pos_chisigma[50]/D");
00270     momInfor->Branch("pos_chisq",pos_chisq,"pos_chisq[50]/D");
00271     momInfor->Branch("neg_chientryNo",neg_chientryNo,"neg_chientryNo[50]/D");
00272     momInfor->Branch("neg_chimean",neg_chimean,"neg_chimean[50]/D");
00273     momInfor->Branch("neg_chimeanerr",neg_chimeanerr,"neg_chimeanerr[50]/D");
00274     momInfor->Branch("neg_chisigma",neg_chisigma,"neg_chisigma[50]/D");
00275     momInfor->Branch("neg_chisq",neg_chisq,"neg_chisq[50]/D");
00276     momInfor->Branch("fitentryNo",fitentryNo,"fitentryNo[50]/D");
00277     momInfor->Branch("fitmean",fitmean,"fitmean[50]/D");
00278     momInfor->Branch("fitmeanerr",fitmeanerr,"fitmeanerr[50]/D");
00279     momInfor->Branch("fitsigma",fitsigma,"fitsigma[50]/D");
00280     momInfor->Branch("fitchisq",fitchisq,"fitchisq[50]/D");
00281     momInfor->Branch("pos_fitentryNo",pos_fitentryNo,"pos_fitentryNo[50]/D");
00282     momInfor->Branch("pos_fitmean",pos_fitmean,"pos_fitmean[50]/D");
00283     momInfor->Branch("pos_fitmeanerr",pos_fitmeanerr,"pos_fitmeanerr[50]/D");
00284     momInfor->Branch("pos_fitsigma",pos_fitsigma,"pos_fitsigma[50]/D");
00285     momInfor->Branch("pos_fitchisq",pos_fitchisq,"pos_fitchisq[50]/D");
00286     momInfor->Branch("neg_fitentryNo",neg_fitentryNo,"neg_fitentryNo[50]/D");
00287     momInfor->Branch("neg_fitmean",neg_fitmean,"neg_fitmean[50]/D");
00288     momInfor->Branch("neg_fitmeanerr",neg_fitmeanerr,"neg_fitmeanerr[50]/D");
00289     momInfor->Branch("neg_fitsigma",neg_fitsigma,"neg_fitsigma[50]/D");
00290     momInfor->Branch("neg_fitchisq",neg_fitchisq,"neg_fitchisq[50]/D");
00291     momInfor->Branch("hits_mean",hits_mean,"hits_mean[50]/D");
00292     momInfor->Branch("hits_sigma",hits_sigma,"hits_sigma[50]/D");
00293     momInfor->Branch("pBin",pBin,"pBin[50]/D");
00294     momInfor-> Branch("ptrk1",ptrk1,"ptrk1[700000]/D");
00295     momInfor-> Branch("dedx1",dedx1,"dedx1[700000]/D");
00296     momInfor->Fill();
00297     momInfor->Write();
00298 
00299     TCanvas c1("c1", "canvas", 500, 400);
00300     c1.Print((m_rootfile+".ps[").c_str());
00301     momInfor -> Draw("dedx1:ptrk1","dedx1>200 && dedx1<1200");
00302     c1.Print((m_rootfile+".ps").c_str());
00303     for(Int_t i=0;i<nbins;i++)
00304     {
00305         m_chi[i]->Draw();
00306         c1.Print((m_rootfile+".ps").c_str());
00307     }
00308     for(Int_t i=0;i<nbins;i++)
00309     {
00310         m_pos_chi[i]->Draw();
00311         c1.Print((m_rootfile+".ps").c_str());
00312     }
00313     for(Int_t i=0;i<nbins;i++)
00314     {
00315         m_neg_chi[i]->Draw();
00316         c1.Print((m_rootfile+".ps").c_str());
00317     }
00318     for(Int_t i=0;i<nbins;i++)
00319     {
00320         m_dedx[i]->Draw();
00321         c1.Print((m_rootfile+".ps").c_str());
00322     }
00323     for(Int_t i=0;i<nbins;i++)
00324     {
00325         m_pos_dedx[i]->Draw();
00326         c1.Print((m_rootfile+".ps").c_str());
00327     }
00328     for(Int_t i=0;i<nbins;i++)
00329     {
00330         m_neg_dedx[i]->Draw();
00331         c1.Print((m_rootfile+".ps").c_str());
00332     }
00333     c1.Print((m_rootfile+".ps]").c_str()); 
00334     f->Close();
00335 
00336     for(int i=0;i<nbins;i++)
00337     {
00338         delete m_chi[i];
00339         delete m_pos_chi[i];
00340         delete m_neg_chi[i];
00341         delete m_dedx[i];
00342         delete m_pos_dedx[i];
00343         delete m_neg_dedx[i];
00344         delete m_hits[i];
00345     }
00346 }