my_star/StRoot/StHbtMaker/Reader/StHbtFlowPicoReader.cxx

#include "StHbtFlowPicoReader.h"
#include "StThreeVectorD.hh"
#include "StHbtMaker/Infrastructure/StHbtEvent.hh"
#include "StHbtMaker/doc/Make/math_constants.h"
#include "SystemOfUnits.h"
#include <fstream>

ClassImp(StHbtFlowPicoReader)

  const double ETAMIN = -5.; /* cutoff values to avoid nan arithmetics */
  const double ETAMAX = 5.;
///////////////////////////////////////////////////////////////////
StHbtFlowPicoReader::StHbtFlowPicoReader(string list_inp,int nevents)
{
  /*
Because pairs analysis is so much slower than DWT, I have to make finer
parallelization splitting (more jobs). For this reason, I have two different
set of "to do list" files: one for the DWT analysis and the other for pairs.
It makes therefore no sense to assume that to_do_list.inp files (for DWT)
will be recycled by the pairs analysis, and consequently to maintain the 
uniformity of the input format. Here I follow the StHbtMuDstReader format
(.lis).
   */
  cout << "StHbtFlowPicoReader::StHbtFlowPicoReader Initializing the reader ! \n";
  to_do_list_inp = list_inp;
  mReaderStatus = 0;           // means "good"
  Nbytes = 0;
  crrnt_entry = 0;
  crrnt_eventI = 0;
  NEvt_in_chain = 0;
  NEVENTS = nevents;
  // read the list of input data files (DST ntuples)


   std::ifstream input(to_do_list_inp.c_str());

   while (input.good())
     { 
       string filename;
       input >>  filename;

       if (filename!="" )
         {
           L_runs.push_back(filename);
         }
     }

   input.close();
cout << " chain has "<<L_runs.size()<<" runs \n";

  Init();
}
///////////////////////////////////////////////////////////////////
int StHbtFlowPicoReader::Init()
{
  cout << " StHbtFlowPicoReader::Init\n"; // open DST

   // build the chain
DST_Chain = new TChain("FlowTree");

   vector<string>::const_iterator it = L_runs.begin();

 
   while (it!=L_runs.end())
     {
       if (NEvt_in_chain <NEVENTS)
         {
           cout << "StHbtFlowPicoReader::StHbtFlowPicoReader add " <<*it<<"\n";

           TChain* tmpChain = new TChain("FlowTree");

           int ok1 = tmpChain->Add((*it).c_str());
           if (!ok1) {cout << "tmpChain->Add error !\n";}
           int tmpEntries = (int)tmpChain->GetEntries();
           NEvt_in_chain += tmpEntries;
           cout << " in file " <<tmpEntries<<" total "<<NEvt_in_chain<< "\n";
           delete tmpChain;
           cout << " adding... \n";
           int ok2 = DST_Chain->Add((*it).c_str());
           if (!ok2) {cout << "DST_Chain->Add error !\n";}
         }
           ++it;
     }

      this->pDST = new flow_pDST(DST_Chain);
       cout << "ST_txtr::ST_txtr flow_pDST instantiated \n";
 this->pDST->fChain->SetBranchStatus("*",0);  // disable all branches

 // enable the branches we need 
 this->pDST->fChain->SetBranchStatus("mRunID",1);
 this->pDST->fChain->SetBranchStatus("mMagneticField",1);
 this->pDST->fChain->SetBranchStatus("mEventID",1);
 this->pDST->fChain->SetBranchStatus("mNtrack",1);
 this->pDST->fChain->SetBranchStatus("mVertexX",1);
 this->pDST->fChain->SetBranchStatus("mVertexY",1);
 this->pDST->fChain->SetBranchStatus("mVertexZ",1);
 this->pDST->fChain->SetBranchStatus("fTracks.fUniqueID",1);
 this->pDST->fChain->SetBranchStatus("fTracks.fBits",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mPtGlobal",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mEtaGlobal",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mPhiGlobal",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mNhits",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mCharge",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mFitPts",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mMaxPts",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mDedx",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mPidElectron",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mPidPion",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mPidKaon",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mPidProton",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mMostLikelihoodPID",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mMostLikelihoodProb",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mDcaGlobal",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mDcaGlobalX",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mDcaGlobalY",1);
 this->pDST->fChain->SetBranchStatus("fTracks.mDcaGlobalZ",1);

 Nentries = (int)(this->pDST->fChain->GetEntries());

 if (Nentries) {return 0;}
 else 
   {
     cout << "StHbtFlowPicoReader::Init() -- see no entries \n";
     return 1;
   }
}
///////////////////////////////////////////////////////////////////
StHbtEvent* StHbtFlowPicoReader::ReturnHbtEvent()
{

  /* using StHbtMcEventReader.cxx as an example
http://www.star.bnl.gov/webdatanfs/dox/html/StHbtMcEventReader_8cxx-source.html

StHbtTrack is here:
http://www.star.bnl.gov/webdatanfs/dox/html/StHbtTrack_8hh-source.html
http://www.star.bnl.gov/webdatanfs/dox/html/StHbtTrack_8cc-source.html

StHbtEvent:
http://www.star.bnl.gov/webdatanfs/dox/html/StHbtEvent_8hh-source.html
http://www.star.bnl.gov/webdatanfs/dox/html/StHbtEvent_8cc-source.html
   */
StHbtEvent* hbtEvent = new StHbtEvent();

 cout << "StHbtFlowPicoReader::ReturnHbtEvent() crrnt_entry begins "
<<crrnt_entry<<"\n";
 int ientry;
 if (crrnt_entry<Nentries)
   {
     ientry = this->pDST->LoadTree(crrnt_entry); 
     nb = this->pDST->fChain->GetEntry(crrnt_entry); Nbytes += nb; 
     
     hbtEvent->SetEventNumber(crrnt_eventI);
     hbtEvent->SetCtbMult(0);
     hbtEvent->SetZdcAdcEast(0);
     hbtEvent->SetZdcAdcWest(0);
     hbtEvent->SetNumberOfTpcHits(0);
     
     int Mult = this->pDST->fTracks_;
     hbtEvent->SetNumberOfTracks(Mult);
     hbtEvent->SetNumberOfGoodTracks(Mult);  
     hbtEvent->SetReactionPlane(0.);
     hbtEvent->SetReactionPlaneSubEventDifference(0.);
     
     float BTesla = pDST->mMagneticField;     
     float vX = pDST->mVertexX;
     float vY = pDST->mVertexY;
     float vZ = pDST->mVertexZ;
     
     StHbtThreeVector VertexPosition = StHbtThreeVector(vX,vY,vZ);

     hbtEvent->SetPrimVertPos(VertexPosition); 
     // loop over tracks

     int UncorrPositivePrim = 0; /* defined in StEventUtilities/StuRefMult.hh*/
     int UncorrNegativePrim = 0;

     for (short itrack = 0; itrack<Mult; itrack++)
       {
	 StHbtThreeVector tmpP;
	 float Phi = pDST->fTracks_mPhiGlobal[itrack];
	 float pT = pDST->fTracks_mPtGlobal[itrack];
	 float eta = pDST->fTracks_mEtaGlobal[itrack];
	 float pX = pT*cos(Phi);
	 float pY = pT*sin(Phi);
	 float expeta = exp(-eta);
	 float theta = 2*atan(expeta);

         short charge = pDST->fTracks_mCharge[itrack];
	 bool reject = false;
	 if (charge ==0) reject = true;
	 /* Sometimes eta is nan. Can we reject such tracks ? */
	 if ( !(eta>ETAMIN) || !(eta<ETAMAX)) reject = true;


	 if (!reject)
	   {
	 StHbtTrack* hbtTrack = new StHbtTrack();
	 hbtTrack->SetHiddenInfo(0);
	     //	     cout << " good eta "<<eta<<"\n";

	     /* zenith angle theta has to be non-negative, between 0 and Pi;
		atan changes between -Pi/2 and Pi/2.
	     */
	     
	     if (theta<0) {theta = C_PI+theta;}
	     float p = pT/sin(theta);
	     float pZ = p-pT*expeta; /* because expeta = tan(theta/2) 
					= (1-cos(theta))/sin(theta) = (p-pZ)/pT */
	     
	     tmpP.setX(pX);
	     tmpP.setY(pY);
	     tmpP.setZ(pZ);
	     hbtTrack->SetP(tmpP); 
	     
	     hbtTrack->SetTrackId(itrack); // ? Frank knows what it is
	     //         int PID=pDST->fTracks_mMostLikelihoodPID[itrack]);
	     float dEdx = pDST->fTracks_mDedx[itrack];

	     float ProbElectron = pDST->fTracks_mPidElectron[itrack];
	     float ProbPion = pDST->fTracks_mPidPion[itrack];
	     float ProbKaon = pDST->fTracks_mPidKaon[itrack];
	     float ProbProton = pDST->fTracks_mPidProton[itrack];
	     
	     hbtTrack->SetdEdx(dEdx);

	     hbtTrack->SetPidProbElectron(ProbElectron);

	     hbtTrack->SetPidProbPion(ProbPion);
	     hbtTrack->SetPidProbKaon(ProbKaon);
	     hbtTrack->SetPidProbProton(ProbProton);


	     int Nhits = pDST->fTracks_mNhits[itrack];
	     
	     float DCA = pDST->fTracks_mDcaGlobal[itrack];
	     float DCAX =  pDST->fTracks_mDcaGlobalX[itrack];
	     float DCAY =  pDST->fTracks_mDcaGlobalY[itrack];
	     float DCAZ =  pDST->fTracks_mDcaGlobalZ[itrack];
	     
	     float DCAXY = sqrt(DCAX*DCAX+DCAY*DCAY);
	     
	     hbtTrack->SetDCAxyGlobal(DCAXY);
	     hbtTrack->SetDCAzGlobal(DCAZ);
	     
	     if (Nhits>=10 && fabs(eta)<0.5 && DCA<3.)
	       {
		 if (charge>0) {UncorrPositivePrim++;}
		 if (charge<0) {UncorrNegativePrim++;}
	       }
	     hbtTrack->SetNHits(charge * Nhits); 
	     //hbtTrack->SetCharge(charge);
	     hbtTrack->SetNHitsPossible(pDST->fTracks_mMaxPts[itrack]);     
	     //hbtTrack->SetPt(pT);
	     //hbtTrack->SetPtGlobal(pT);

	     int gid =  pDST->fTracks_mMostLikelihoodPID[itrack];
	     float pid_prob = pDST->fTracks_mMostLikelihoodProb[itrack];

	     // Assume that this can only happen in MC
	     if (pid_prob == 1.)
	       {
		 if (gid ==2 || gid ==3) 
		   {
		     hbtTrack->SetNSigmaElectron(0.); 
		     hbtTrack->SetNSigmaPion(-1000.); 
		     hbtTrack->SetNSigmaKaon(-1000.); 
		     hbtTrack->SetNSigmaProton(-1000.); 
		   }
		 else if (gid==8 || gid ==9)
		   {
		     hbtTrack->SetNSigmaElectron(1000.); 
		     hbtTrack->SetNSigmaPion(0.); 
		     hbtTrack->SetNSigmaKaon(-1000.); 
		     hbtTrack->SetNSigmaProton(-1000.); 
		   }
		 else if (gid==11 || gid==12)
		   {
		     hbtTrack->SetNSigmaElectron(1000.); 
		     hbtTrack->SetNSigmaPion(1000.); 
		     hbtTrack->SetNSigmaKaon(0.); 
		     hbtTrack->SetNSigmaProton(-1000.); 
		   }
		 else if (gid==14 ||  gid==15)
		   {
		     hbtTrack->SetNSigmaElectron(1000.); 
		     hbtTrack->SetNSigmaPion(1000.); 
		     hbtTrack->SetNSigmaKaon(1000.); 
		     hbtTrack->SetNSigmaProton(0.); 
		   }
		 else
		   {
		     hbtTrack->SetNSigmaElectron(1000.); 
		     hbtTrack->SetNSigmaPion(0.); 
		     hbtTrack->SetNSigmaKaon(-1000.); 
		     hbtTrack->SetNSigmaProton(-1000.); 
		   }

	       }
	     
	     // define helix:
	     
	     const StPhysicalHelixD& mHelix = 
	       StPhysicalHelixD(tmpP,VertexPosition,BTesla*tesla,charge);
	     
	     hbtTrack->SetHelix(mHelix);
	     
	     hbtEvent->TrackCollection()->push_back(hbtTrack);
	   }
       }

     hbtEvent->SetUncorrectedNumberOfPositivePrimaries(UncorrPositivePrim);
     hbtEvent->SetUncorrectedNumberOfNegativePrimaries(UncorrNegativePrim);

     crrnt_entry++;
     crrnt_eventI++;
          cout << "StHbtFlowPicoReader return event # "<<crrnt_eventI<<"\n";
	  //	  hbtEvent->RandomizeTrackOrder();
     return hbtEvent;
   }
 else
   {
     mReaderStatus = 1;
     delete hbtEvent;
     return 0;
   }
}
///////////////////////////////////////////////////////////////////
/*
int StHbtFlowPicoReader::WriteHbtEvent(StHbtEvent*)
{
  cout << " calling my own  StHbtFlowPicoReader::WriteHbtEvent(StHbtEvent*)\n";
  return 0;
}
*/
///////////////////////////////////////////////////////////////////

StHbtString StHbtFlowPicoReader::Report()
{
return (StHbtString)" Here is the reader report \n";
}
///////////////////////////////////////////////////////////////////
StHbtFlowPicoReader::~StHbtFlowPicoReader()
{
}