StRoot: StHighptPool/StHiMicroAnalysis/StHiSpectra.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: StHiSpectra.cxx,v 1.4 2003/04/30 20:37:30 perev Exp $                                    
00004  *
00005  * Author: Bum Choi, UT Austin, Apr 2002
00006  *
00007  ***************************************************************************
00008  *
00009  * Description:  Class for making highpt inclusive spectra from highpt
00010  *               uDST's
00011  *
00012  ***************************************************************************
00013  * 
00014  * $Log: StHiSpectra.cxx,v $
00015  * Revision 1.4  2003/04/30 20:37:30  perev
00016  * Warnings cleanup. Modified lines marked VP
00017  *
00018  * Revision 1.3  2002/05/31 21:58:29  jklay
00019  * Updated analysis code to use new cut class
00020  *
00021  * Revision 1.2  2002/04/03 00:23:27  jklay
00022  * Fixed private member access bugs in analysis code
00023  *
00024  * Revision 1.1  2002/04/02 20:05:18  jklay
00025  * Bums analysis tools for highpt uDSTs
00026  *
00027  *
00028  **************************************************************************/
00029 #include "StHiSpectra.h"
00030 
00031 
00032 const char* border = "*************************************************";
00033 
00034 //__________________
00035 
00036 StHiSpectra::StHiSpectra(const char* inputDir,
00037                          const char* outRootName,
00038                          const char* f)
00039   : StHiBaseAnalysis(inputDir,outRootName)
00040      
00041 {
00042 
00043   mEfficiencyString = ""; mEfficiencyFileName = f;
00044   mEfficiencyMap = 0;
00045 }
00046 
00047 //__________________
00048 
00049 StHiSpectra::~StHiSpectra()
00050 {
00051 }
00052 
00053 //__________________
00054 //
00055 // reads in all the files in the directory
00056 // by default, looks for files ending in minimc.root 
00057 //
00058 
00059 Int_t
00060 StHiSpectra::initMore()
00061 {
00062 
00063   //  char name[100];
00064   Int_t stat=0;
00065 
00066   TString fName = "fEfficiency";
00067   TFile* file=0; 
00068   if(mEfficiencyFileName==""){
00069     cout << "no efficiency file set " << endl;
00070   }
00071   else{
00072     file=new TFile(mEfficiencyFileName.Data());
00073     if(!file ||!file->IsOpen()){
00074       cout << "Cannot find efficiency file : " << mEfficiencyFileName << endl;
00075     }
00076     else{
00077         mEfficiencyMap = (TF2*)file->Get(fName.Data());
00078         delete file;
00079     }
00080   }
00081   if(!mEfficiencyMap){
00082     cout << "Cannot find efficency function : " << fName.Data() <<endl;
00083     cout << "WARNING-USING DUMMY EFFIENCY" << endl;
00084     mEfficiencyMap = new TF2("fEfficiency","1");
00085   }
00086   
00087   
00088   cout << "\tefficiency file : " << mEfficiencyFileName << endl;
00089   cout << "\tefficiency fcn title : " << mEfficiencyMap->GetTitle() << endl;
00090 
00091   /*
00092   if(mEfficiencyString==""){
00093     cout << "no efficiency string" << endl;
00094     stat++;
00095   }
00096   
00097   else{
00098     // create the efficiency function
00099     // x is eta, y is pt
00100     cout << "efficiency : " << mEfficiencyString << endl;
00101     
00102     mEfficiencyMap = new TF2("fEfficiencyMap",mEfficiencyString,-.55,.55,1.6,6);
00103   }
00104   */
00105   return stat;
00106 }
00107 
00108 
00109 void
00110 StHiSpectra::initHistograms()
00111 {
00112   cout << "StHiSpectra::initHistograms()" << endl;
00113 
00114   //*********************
00115 
00116   using namespace Bin;
00117 
00118   //*********************
00119 
00120   //  gStyle->SetPalette(1,0);
00121 
00122   char name[500];
00123 
00124   h1CentralityCut = new TH1D("h1CentralityCut","centrality cut",
00125                              nFlowCentBin,flowCentMin,flowCentMax);
00126 
00127   h1Centrality = new TH1D("h1Centrality","centrality",
00128                              nFlowCentBin,flowCentMin,flowCentMax);
00129 
00130   h1VertexZ = new TH1D("h1VertexZ","vertex z",
00131                        nVertexZEvtBin,vertexZEvtMin,vertexZEvtMax);
00132 
00133 //  sprintf(title,"vertex z (%d<=zdc cent<=%d)",
00134 //        CutRc::mZdcCtbCent[1],CutRc::mZdcCtbCent[0]);
00135 //  h1VertexZCut = new TH1D("h1VertexZCut",title,
00136 //                        nVertexZEvtBin,vertexZEvtMin,vertexZEvtMax);
00137 
00138   h1VertexZThin = new TH1D("h1VertexZThin","vertex z",
00139                            nVertexZEvtThinBin,vertexZEvtMin,vertexZEvtMax);
00140 
00141 //  sprintf(title,"vertex r %d<=zdc cent<=%d)",
00142 //        CutRc::mZdcCtbCent[1],CutRc::mZdcCtbCent[0]);
00143 //  h1VertexRCut = new TH1D("h1VertexRCut",title,
00144 //                        nVertexREvtBin,vertexREvtMin,vertexREvtMax);
00145 
00146   cout << "using bins" << endl;
00147   for(int iBin=0; iBin<2; iBin++){
00148     varBin[iBin].ptBinsAry = new TArrayD;
00149     initPtAry(varBin[iBin].ptBinsAry,iBin);
00150     int nBin = varBin[iBin].ptBinsAry->GetSize()-1;
00151     varBin[iBin].nPtBinAry = nBin;;
00152     cout << "Bin " << iBin << " : " << endl;
00153     for(int i=0; i<nBin; i++){
00154       cout << varBin[iBin].ptBinsAry->At(i) << ", ";
00155     }
00156     cout << varBin[iBin].ptBinsAry->At(nBin) << endl;
00157   }
00158   //------------------------------------------------------------
00159 
00160   TString sPM[2] = { "Plus","Minus"};
00161 
00162   //********************************************************
00163   // number of events used to scale the spectra
00164 
00165   h1NEvent = new TH1D("h1NEvent","h1NEvent",1,1,2);
00166 
00167   // eta cut
00168 
00169   h1EtaCut = new TH1D("h1EtaCut","h1EtaCut",2,0,2);
00170 
00171 
00172   //********************************************************
00173   // yield for both signs
00174 
00175   h2Yield
00176     = new TH2D("h2Yield","h2Yield",
00177                nEtaSmallBin,etaSmallMin,etaSmallMax,
00178                nPtBin,ptMin,ptMax);
00179   h2Yield->Sumw2();
00180 
00181   
00182   
00183 
00184   //********************************************************
00185   // var bin
00186   
00187   for(Int_t iBin=0; iBin<mNVarBin; iBin++){
00188 
00189     //**** both charge signs
00190 
00191     setName(name,"h1OneOverPt",iBin);
00192     varBin[iBin].mean.h1OneOverPt =
00193       new TH1D(name,name,varBin[iBin].nPtBinAry,
00194                varBin[iBin].ptBinsAry->GetArray());
00195     varBin[iBin].mean.h1OneOverPt->Sumw2();
00196     
00197     setName(name,"h1WeightedMean",iBin);
00198     varBin[iBin].mean.h1WeightedMean =
00199       new TH1D(name,name,varBin[iBin].nPtBinAry,
00200                varBin[iBin].ptBinsAry->GetArray());
00201     varBin[iBin].mean.h1WeightedMean->Sumw2();
00202 
00203 
00204     //
00205     // raw
00206     //
00207     setName(name,"h1Raw",iBin);
00208     varBin[iBin].spec.h1Raw =
00209       new TH1D(name,name,varBin[iBin].nPtBinAry,
00210                varBin[iBin].ptBinsAry->GetArray());
00211     varBin[iBin].spec.h1Raw->Sumw2();
00212     
00213     //
00214     // efficiency corrected
00215     //
00216     setName(name,"h1EffCorrected",iBin);
00217     varBin[iBin].spec.h1EffCorrected =
00218       new TH1D(name,name,varBin[iBin].nPtBinAry,
00219                varBin[iBin].ptBinsAry->GetArray());
00220     varBin[iBin].spec.h1EffCorrected->Sumw2();
00221     
00222     //
00223     // corrected (background + momentum resolution)
00224     //
00225     setName(name,"h1Corrected",iBin);
00226     varBin[iBin].spec.h1Corrected =
00227       new TH1D(name,name,varBin[iBin].nPtBinAry,
00228                varBin[iBin].ptBinsAry->GetArray());
00229     varBin[iBin].spec.h1Corrected->Sumw2();
00230 
00231 
00232     //**** spectra plus, minus
00233 
00234     for(Int_t iCharge=0; iCharge<2; iCharge++){
00235 
00236       // mean stuff
00237       //
00238 
00239       setName(name,"h1OneOverPt",iBin,sPM[iCharge].Data());
00240       varBin[iBin].meanPM[iCharge].h1OneOverPt =
00241         new TH1D(name,name,varBin[iBin].nPtBinAry,
00242                  varBin[iBin].ptBinsAry->GetArray());
00243       varBin[iBin].meanPM[iCharge].h1OneOverPt->Sumw2();
00244 
00245       setName(name,"h1WeightedMean",iBin,sPM[iCharge].Data());
00246       varBin[iBin].meanPM[iCharge].h1WeightedMean =
00247         new TH1D(name,name,varBin[iBin].nPtBinAry,
00248                  varBin[iBin].ptBinsAry->GetArray());
00249       varBin[iBin].meanPM[iCharge].h1WeightedMean->Sumw2();
00250 
00251       //
00252       // raw
00253       //
00254       setName(name,"h1Raw",iBin,sPM[iCharge].Data());
00255       varBin[iBin].specPM[iCharge].h1Raw =
00256         new TH1D(name,name,varBin[iBin].nPtBinAry,
00257                  varBin[iBin].ptBinsAry->GetArray());
00258       varBin[iBin].specPM[iCharge].h1Raw->Sumw2();
00259 
00260       //
00261       // efficiency corrected
00262       //
00263       setName(name,"h1EffCorrected",iBin,sPM[iCharge].Data());
00264       varBin[iBin].specPM[iCharge].h1EffCorrected =
00265         new TH1D(name,name,varBin[iBin].nPtBinAry,
00266                  varBin[iBin].ptBinsAry->GetArray());
00267       varBin[iBin].specPM[iCharge].h1EffCorrected->Sumw2();
00268       
00269       //
00270       // corrected (background + momentum resolution)
00271       //
00272       
00273       setName(name,"h1Corrected",iBin,sPM[iCharge].Data());
00274       varBin[iBin].specPM[iCharge].h1Corrected =
00275         new TH1D(name,name,varBin[iBin].nPtBinAry,
00276                  varBin[iBin].ptBinsAry->GetArray());
00277       varBin[iBin].specPM[iCharge].h1Corrected->Sumw2();
00278       
00279       
00280     }
00281 
00282     //***** ratios
00283 
00284     setName(name,"h1RawRatio",iBin);
00285     varBin[iBin].h1RawRatio
00286       = new TH1D(name,name,varBin[iBin].nPtBinAry,
00287                  varBin[iBin].ptBinsAry->GetArray());
00288     varBin[iBin].h1RawRatio->Sumw2();
00289 
00290     setName(name,"h1CorrectedRatio",iBin);
00291     varBin[iBin].h1CorrectedRatio
00292       = new TH1D(name,name,varBin[iBin].nPtBinAry,
00293                  varBin[iBin].ptBinsAry->GetArray());
00294     varBin[iBin].h1CorrectedRatio->Sumw2();
00295 
00296   }  // varBin
00297 }
00298 
00299 //_____________________
00300 
00301 void 
00302 StHiSpectra::trackLoop()
00303 {
00304   if(mDebug)
00305     cout << "StHiSpectra::spectraTrackLoop()" << endl;
00306 
00307   Int_t nTrack = mHiMicroEvent->NTrack();
00308   StHiMicroTrack* track;
00309   
00310 //VPunused  Float_t vertexZ = mHiMicroEvent->VertexZ();
00311 
00312   for(Int_t i=0; i<nTrack; i++){
00313     track =(StHiMicroTrack*) mHiMicroEvent->tracks()->At(i);
00314 
00315     //***************** spectra  ***********************************
00316     
00317     if(!CutRc::AcceptTrackHalf(track)) {
00318         continue; }
00319     if(!CutRc::Accept(track)) continue; 
00320 
00321     //**************************************************************
00322     Float_t ptPr = track->PtPr();
00323     Float_t eta = track->EtaPr();
00324 //VPunused    Float_t ptGl = track->PtGl();
00325     //Float_t midZ = 100*TMath::Tan(track->DipAnglePr()) + vertexZ ;
00326     //Float_t exitZ = 200*TMath::Tan(track->DipAnglePr()) + vertexZ ;
00327     //
00328     // efficiency correction here
00329     //
00330     Float_t eff = mEfficiencyMap->Eval(eta,ptPr);
00331     
00332     Float_t weight = (1./eff);
00333 
00334     //
00335     // yield for both charge signs
00336 
00337     h2Yield->Fill(eta,ptPr);
00338 
00339     //
00340     // counts
00341     //
00342     if(ptPr>1.5&&ptPr<6){
00343       mNHiPtTrack++;
00344       if(mNHiPtTrack%10000==0){
00345         cout << "pt: " << ptPr << " eta: " << eta << " eff : " << eff << endl;
00346       }
00347     }
00348     //
00349     // spectra stuff here
00350     //
00351 
00352     Int_t iCharge = (track->Charge()>0) ? 0 : 1; //plus is 0
00353 
00354     for(Int_t iBin=0; iBin<mNVarBin; iBin++){
00355       
00356       //
00357       // both
00358       //
00359       
00360       varBin[iBin].spec.h1Raw->Fill(ptPr);
00361       varBin[iBin].mean.h1OneOverPt->Fill(ptPr,1./ptPr);
00362       varBin[iBin].spec.h1EffCorrected->Fill(ptPr,weight);
00363       varBin[iBin].spec.h1Corrected->Fill(ptPr,weight);
00364      
00365 
00366       //
00367       // plus/minus
00368       //
00369 
00370       // mean stuff
00371       //
00372       varBin[iBin].meanPM[iCharge].h1OneOverPt->Fill(ptPr,1./ptPr);
00373 
00374       varBin[iBin].specPM[iCharge].h1Raw->Fill(ptPr);
00375       varBin[iBin].specPM[iCharge].h1EffCorrected->Fill(ptPr,weight);
00376       varBin[iBin].specPM[iCharge].h1Corrected->Fill(ptPr,weight);
00377 
00378       // vertex z east, west - make sure they stay on the same side
00379       //
00380 
00381     } // varbin
00382 
00383   }
00384 }
00385 
00386 //_____________________
00387 
00388 void
00389 StHiSpectra::fillEventHistograms()
00390 {
00391   h1Centrality->Fill(mHiMicroEvent->Centrality());
00392 
00393   h1VertexZ->Fill(mHiMicroEvent->VertexZ());
00394                   
00395   h1VertexZThin->Fill(mHiMicroEvent->VertexZ());
00396 
00397   if(CutRc::AcceptVertexZ(mHiMicroEvent))
00398     h1CentralityCut->Fill(mHiMicroEvent->Centrality());
00399 
00400 //  if(CutRc::AcceptZdcCtbCent(mHiMicroEvent)){
00401 //    h1VertexZCut->Fill(mHiMicroEvent->VertexZ());
00402 //    //h1VertexRCut->Fill(mHiMicroEvent->VertexR());
00403 //  }
00404 }
00405 
00406 
00407 //______________________
00408 
00409 Bool_t
00410 StHiSpectra::acceptEvent(StHiMicroEvent* event)
00411 {
00412   return CutRc::Accept(event);
00413 
00414 }
00415 
00416 
00417 //______________________
00418 
00419 void
00420 StHiSpectra::finishHistograms()
00421 {
00422   h1NEvent->SetBinContent(1,mNEventAccepted);
00423   h1EtaCut->SetBinContent(1,CutRc::mEta[0]);
00424   h1EtaCut->SetBinContent(2,CutRc::mEta[1]);
00425 
00426   for(Int_t iBin=0; iBin<mNVarBin; iBin++){
00427 
00428     varBin[iBin].mean.h1WeightedMean->Divide(varBin[iBin].spec.h1Raw,varBin[iBin].mean.h1OneOverPt);
00429     
00430     for(Int_t iCharge=0; iCharge<2; iCharge++){
00431 
00432       varBin[iBin].meanPM[iCharge].h1WeightedMean->Divide(varBin[iBin].specPM[iCharge].h1Raw,varBin[iBin].meanPM[iCharge].h1OneOverPt);
00433 
00434     }
00435   }
00436 
00437 }
00438 
00439 
00440 
00441 //_____________________________________
00442 
00443 ClassImp(StHiSpectra)