Nikita Ermakov 2 years ago
parent
commit
ed00932b7c
35 changed files with 7539 additions and 11 deletions
  1. 3 3
      StRoot/StFemtoDstMaker/StFemtoDstPythia6Maker.cxx
  2. 461 0
      StRoot/StFemtoDstMaker/StHbtStarGenHijingDstReader.cxx
  3. 77 0
      StRoot/StFemtoDstMaker/StHbtStarGenHijingDstReader.h
  4. 451 0
      macros/HBT/expPionDCHBT_pau200_y2015_3D.C
  5. 444 0
      macros/HBT/simPionDCHBT_pau200_3D.C
  6. 194 0
      macros/femtoDst/FemtoDstMaker_mesons_pp500_y2011.C
  7. 1 1
      macros/femtoDst/FemtoDstMaker_pions_pp500_y2011.C
  8. 63 0
      macros/femtoDst/FemtoDstQA_pp500_y2011.C
  9. 6 7
      macros/simulation/starHijing.C
  10. 157 0
      source/HistoCollector/BuildDependencies.C
  11. 57 0
      source/HistoCollector/ComparePID.C
  12. 575 0
      source/HistoCollector/DataProcessing.cxx
  13. 140 0
      source/HistoCollector/DataProcessing.h
  14. 111 0
      source/HistoCollector/FemtoUtils.h
  15. 1643 0
      source/HistoCollector/HbtFitter.cxx
  16. 496 0
      source/HistoCollector/HbtFitter.h
  17. 297 0
      source/HistoCollector/HistoCollector1D.cxx
  18. 64 0
      source/HistoCollector/HistoCollector1D.h
  19. 210 0
      source/HistoCollector/HistoCollector3D.cxx
  20. 50 0
      source/HistoCollector/HistoCollector3D.h
  21. 16 0
      source/HistoCollector/PlotFinalResults.C
  22. 86 0
      source/HistoCollector/PresentAllEnergies.C
  23. 924 0
      source/HistoCollector/PresentResults.C
  24. 30 0
      source/HistoCollector/ProcessMesons1D.C
  25. 624 0
      source/HistoCollector/ProcessMesons3D.C
  26. 51 0
      source/HistoCollector/Processor.C
  27. 27 0
      source/HistoCollector/RunProcessor.C
  28. 120 0
      source/HistoCollector/Test3DFitter.C
  29. 24 0
      xml/HBT/expPionDCHBT_pau200_y2015_3D.xml
  30. 26 0
      xml/femtoDst/FemtoDstMaker_mesons_pp500_y2011.xml
  31. 26 0
      xml/femtoDst/FemtoDstMaker_mesons_pp500_y2011_p1.xml
  32. 26 0
      xml/femtoDst/FemtoDstMaker_mesons_pp500_y2011_p2.xml
  33. 26 0
      xml/femtoDst/FemtoDstMaker_mesons_pp500_y2011_p3.xml
  34. 21 0
      xml/femtoDst/FemtoDstQA_pp500_y2011.xml
  35. 12 0
      xml/others/sum.xml

+ 3 - 3
StRoot/StFemtoDstMaker/StFemtoDstPythia6Maker.cxx

@@ -263,11 +263,11 @@ Int_t StFemtoDstPythia6Maker::Make() {
       //Default set to the first bit: primary-vertex-used
       tmap1 |= 1 << 0;
     }
-    else if(ptot < 0.3) { //Tracks may come to the end of the outer sector: R = 2m
+    else if (ptot < 0.3) { //Tracks may come to the end of the outer sector: R = 2m
       tmap1 = 0;
       tmap2 = 0;
-      for(int iBit=0; iBit<32; iBit++) {
-
+      for (int iBit=0; iBit<32; iBit++)
+      {
 	if(randy->Rndm() < 0.95) {
 	  tmap1 |= 1 << iBit;
 	}

+ 461 - 0
StRoot/StFemtoDstMaker/StHbtStarGenHijingDstReader.cxx

@@ -0,0 +1,461 @@
+#include "StHbtStarGenHijingDstReader.h"
+#include "TString.h"
+
+ClassImp(StHbtStarGenHijingDstReader)
+
+//_________________
+StHbtStarGenHijingDstReader::StHbtStarGenHijingDstReader(const char *dir, 
+							 const char *fileName,
+							 const char *filter,
+							 int maxFiles,
+							 bool debug) {
+    mDir      = string(dir);
+    mFileName = string(fileName);
+    mFilter   = filter;
+    mMaxFiles = maxFiles;
+    mDebug    = debug;
+
+    mTotalTracks = 0.;
+  
+    mTree       = 0;
+    mTChain     = 0;
+    mSGEvent    = 0;
+    mEventIndex = 0;
+
+    randy = new TRandom3(0);
+
+    InitRead(mDir, mFileName, mFilter, mMaxFiles);
+    mNEvents = (unsigned int)mTChain->GetEntries();
+}
+
+//_________________
+int StHbtStarGenHijingDstReader::FillChain(TChain *chain, char *dir,
+					   const char *filter, int maxFiles) {
+
+    TSystem *gSystem = new TSystem();
+    if (!gSystem) {
+	std::cout << "StHbtStarGenHijingDstReader[ERROR]: cannot allocate memory for TSystem\n";
+	return 0;
+    }
+    void *lDir = gSystem->OpenDirectory(dir);
+    if (!lDir) {
+	std::cout << "StHbtStarGenHijingDstReader[ERROR]: can't open directory " 
+		  << dir << std::endl;
+	delete gSystem;
+	return 0;
+    }
+  
+    int mCount = 0;
+    char *lFileName;
+
+    while ((lFileName = (char *)gSystem->GetDirEntry(dir))) {
+	if (strcmp(lFileName, ".") == 0 || strcmp(lFileName, "..") == 0)
+	    continue;
+	if (strstr(lFileName, filter)) {
+	    char *lFullName = gSystem->ConcatFileName(dir, lFileName);
+	    std::cout << "StHbtStarGenHijingDstReader[INFO]: Adding "
+		      << lFullName << " to the chain" << std::endl;
+	    chain->Add(lFullName);
+	    delete lFullName;
+	    mCount++;
+	    if ((maxFiles > 0) && (mCount > maxFiles)) break;
+	}
+    }
+  
+    std::cout << "StHbtStarGenHijingDstReader[INFO]: Added " << mCount 
+	      << " files to the chain" << std::endl;
+    delete gSystem;
+    return mCount;
+}
+
+//_________________
+int StHbtStarGenHijingDstReader::FillChain(TChain *chain, const char *fileName,
+					   int maxFiles) {
+
+    std::ifstream lInStream(fileName);
+
+    if (!lInStream.is_open()) {
+	std::cout << "StHbtStarGenHijingDstReader[ERROR]: can't open file list: "
+		  << fileName << std::endl;
+	return 0;
+    }
+    std::cout << "StHbtStarGenHijingDstReader[INFO]: using file list: "
+	      << fileName << std::endl;
+
+    int mCount = 0;
+    char buf[0xFF];
+    for ( ; lInStream.good(); ) {
+	lInStream.getline(buf, 0xFF);
+	TString lFileName(buf);
+	if (lFileName.Contains("root")) {
+	    chain->Add(buf);
+	    mCount++;
+	    if ((maxFiles > 0) && (mCount > maxFiles)) break;
+	}
+    }
+
+    std::cout << "StHbtStarGenHijingDstReader[INFO]: Added " << mCount
+	      << " files to the chain" << std::endl;
+    return mCount;
+}
+
+//_________________
+int StHbtStarGenHijingDstReader::InitRead(string aDir, string fileName,
+					  string aFilter, int maxFiles) {
+    mEventIndex = 0;
+    //  mSGEvent = new StSGEvent();
+    mTChain = new TChain("StFemtoDst", "StFemtoDst");
+  
+    std::cout << "StHbtStarGenHijingDstReader[INFO]: Call InitRead" << std::endl;
+  
+    int lNumFiles = 0;
+    if (!fileName.empty()) {
+	if (strstr(fileName.c_str(), ".lis") ||
+	    strstr(fileName.c_str(), ".list")) {
+	    lNumFiles = FillChain( mTChain, 
+				   (aDir + fileName).c_str(), 
+				   maxFiles );
+	}
+	else {
+	    mTChain->Add((aDir + fileName).c_str());
+	    lNumFiles++;
+	}
+    } else {
+	lNumFiles = FillChain(mTChain, (char *)aDir.c_str(), // not cool
+			      aFilter.c_str(), maxFiles);
+    }
+    mTChain->SetBranchAddress("StarGenEvent", &mSGEvent);
+    mNEvents = (unsigned int)mTChain->GetEntries();
+    std::cout << "StHbtStarGenHijingDstReader::InitRead [INFO]"
+	      << " NEvents to read: " << mNEvents << std::endl;
+    return lNumFiles;
+}
+
+//_________________
+void StHbtStarGenHijingDstReader::UninitRead() {
+    //  if (mSGEvent) delete mSGEvent;
+    if (mTChain) delete mTChain;
+    mSGEvent = 0;
+    mTChain   = 0;
+}
+
+//_________________
+void StHbtStarGenHijingDstReader::PrintTotalTracks() {
+
+    std::cout << "StHbtStarGenHijingDstReader[INFO]: mTotalTracks = "
+	      << mTotalTracks << std::endl;
+}
+
+//_________________
+int StHbtStarGenHijingDstReader::GetNEvents() {
+
+    if (mTChain) {
+	return mNEvents;
+    }
+    else {
+	return -1;
+    }
+}
+
+//_________________
+StHbtEvent *StHbtStarGenHijingDstReader::Read() {
+
+    if(mDebug) {
+	std::cout << "StHbtStarGenHijingDstReader[DEBUG]: mEventIndex = "
+		  << mEventIndex << " mNEvents = "
+		  << mNEvents << std::endl;
+    }
+      
+    if(!mNEvents) {
+	std::cout << "StHbtStarGenHijingDstReader[INFO]: no events to read\n";
+	return 0;
+    }
+
+    mSGEvent->Clear();
+    int lBytes = mTChain->GetEntry(mEventIndex++);
+
+    if(mNEvents < mEventIndex) {
+	std::cout << "StHbtStarGenHijingDstReader[INFO]: EOF\n";
+	return 0;
+    }
+
+    if(!lBytes) {
+	std::cout << "StHbtStarGenHijingDstReader[INFO]: no event\n";
+	return 0;
+    }
+
+    StHbtEvent *mHbtEvent = NULL;
+    
+    if (mSGEvent)
+    {
+	mHbtEvent = new StHbtEvent(); // create new StHbtEvent
+	const float magRFF = -4.98;   // [10^3 G], ref: p+Au 200 GeV y2015
+	StThreeVectorF v(0., 0., 0.); // event vertex
+	int nTracks = mSGEvent->GetNumberOfParticles(); // number of tracks in an event
+
+	//
+	// Set an event information.
+	//
+	//  Number of the positive and negative tracks
+	//  are setting after the loop over tracks.
+	//
+	mHbtEvent->SetEventNumber(mSGEvent->GetEventNumber());
+	mHbtEvent->SetRunNumber(mSGEvent->GetRunNumber());
+	mHbtEvent->SetZdcAdcEast(0.);
+	mHbtEvent->SetZdcAdcWest(0.);
+	mHbtEvent->SetUncorrectedNumberOfPrimaries(nTracks);
+	mHbtEvent->SetPrimVertPos(v);
+	mHbtEvent->SetMagneticField(magRFF);
+	mHbtEvent->SetVpdVz(0.);
+	mHbtEvent->SetCent9(0);
+	mHbtEvent->SetCent16(0);
+	mHbtEvent->SetRefMultCorr(nTracks);
+	mHbtEvent->SetSimulation(true);
+
+	mTotalTracks += (float)nTracks;
+	if (mDebug)
+	    std::cout << "StHbtStarGenHijingDstReader[DEBUG]: nTracks = "
+		      << nTracks << std::endl;
+
+	int nTracksPos = 0;
+	int nTracksNeg = 0;
+
+	for (int i = 0; i < nTracks; i++)
+	{
+	    StarGenParticle *sParticle = (*mSGEvent)[i];
+	    StHbtTrack *mHbtTrack = new StHbtTrack;
+
+	    //
+	    // Calculate charge of the particle
+	    //
+	    int charge;
+	    int pdg = sParticle->GetId();
+		    
+	    if (pdg >= 0)
+	    {
+		nTracksPos++;
+		charge = +1;
+	    }
+	    else
+	    {
+		nTracksNeg++;
+		charge = -1;
+	    }
+
+	    mHbtTrack->SetCharge(charge);
+	    mHbtTrack->SetTrackId(i);
+	    mHbtTrack->SetTrackType(1);   // 1 - primary track
+	    mHbtTrack->SetFlag(500);      // old flag < 1000 means no pile up
+
+	    //
+	    // Energy loss of the particle
+	    //
+	    int nHits = 45;
+	    const TLorentzVector &p = sParticle->momentum();
+		    
+	    mHbtTrack->SetNHits(nHits*charge);
+	    mHbtTrack->SetNHitsFit(nHits);
+	    mHbtTrack->SetNHitsPossible(nHits);
+	    mHbtTrack->SetNHitsDedx(nHits);
+	    mHbtTrack->SetdEdx(dedxMean(sParticle->GetMass(), p.P()));
+	    mHbtTrack->SetChiSquaredXY(1.);
+	    mHbtTrack->SetChiSquaredZ(1.);
+	    
+	    //
+	    // Calculate nSigma of the
+	    //  electron, pion, kaon, proton
+	    //
+	    float nSigmaElectron = -999.;
+	    float nSigmaPion = -999.;
+	    float nSigmaKaon = -999.;
+	    float nSigmaProton = -999.;
+	    
+	    switch (abs(pdg))
+	    {
+	    case 211:
+		nSigmaPion = 0.;
+		break;
+	    case 321:
+		nSigmaKaon = 0.;
+		break;
+	    case 2212:
+		nSigmaProton = 0.;
+		break;
+	    }
+	    
+	    mHbtTrack->SetNSigmaElectron(nSigmaElectron);
+	    mHbtTrack->SetNSigmaPion(nSigmaPion);
+	    mHbtTrack->SetNSigmaKaon(nSigmaKaon);
+	    mHbtTrack->SetNSigmaProton(nSigmaProton);
+
+	    //
+	    // Calculate distance of closet approach
+	    //
+	    float xDCA = sParticle->GetVx();
+	    float yDCA = sParticle->GetVy();
+	    float zDCA = sParticle->GetVz();
+	    float xyDCA = sqrt(xDCA*xDCA + yDCA*yDCA);
+	    
+	    mHbtTrack->SetDCAxy(xyDCA);
+	    mHbtTrack->SetDCAz(zDCA);
+	    mHbtTrack->SetDCAxyGlobal(xyDCA);
+	    mHbtTrack->SetDCAzGlobal(zDCA);
+
+	    //
+	    // Set momentum of the particle
+	    //
+	    StThreeVectorF mGlobMomentum(p.Px(),
+					 p.Py(),
+					 p.Pz());
+	    StThreeVectorF mPrimMomentum(mGlobMomentum);
+	    mHbtTrack->SetP(mPrimMomentum);
+
+	    if(mDebug)
+		std::cout << Form("Track parameters: px=%3.2f py=%3.2f pz=%3.2f q=%1d nhits=%2d",
+				  mPrimMomentum.x(), mPrimMomentum.y(), mPrimMomentum.z(),
+				  charge, nHits) << std::endl;
+	    
+	    //
+	    // Calculate helix
+	    //
+	    StThreeVectorF mGlobOrigin(xDCA, yDCA, zDCA);
+	    StThreeVectorF mPrimOrigin(0., 0., 0.);
+	    StPhysicalHelixD mGlobHelix(mGlobMomentum, mGlobOrigin, magRFF*kilogauss, charge);
+	    StPhysicalHelixD mPrimHelix(mPrimMomentum, mPrimOrigin, magRFF*kilogauss, charge);
+
+	    mHbtTrack->SetHelixGlobal(mGlobHelix);
+	    mHbtTrack->SetHelix(mPrimHelix);
+
+	    //
+	    // Calculate topology map
+	    //
+	    unsigned int tmap1 = 0;
+	    unsigned int tmap2 = 0;
+	    float ptot = p.P();
+	    if (ptot < 0.15) // tracks only in the inner sector. R=1m
+	    {
+		tmap1 = 0;
+		tmap2 = 0;
+		for (int iBit = 0; iBit < 32; iBit++)
+		{
+		    if (randy->Rndm() < 0.95)
+			tmap1 |= 1 << iBit;
+		    else
+			tmap1 |= 0 << iBit;
+		}
+		tmap1 |= 1 << 0; // default set to the first bit: primary-vertex-used
+	    }
+	    else
+	    {
+		if (ptot < 0.3) // tracks may come to the end of the outer sector: R = 2m
+		{ 
+		    tmap1 = 0;
+		    tmap2 = 0;
+		    for (int iBit = 0; iBit < 32; iBit++)
+		    {
+			if (randy->Rndm() < 0.95)
+			    tmap1 |= 1 << iBit;
+			else 
+			    tmap1 |= 0 << iBit;
+			
+			if (iBit < 27) 
+			    if (randy->Rndm() < 0.85) 
+				tmap2 |= 1 << iBit;
+			    else 
+				tmap2 |= 0 << iBit;
+		    }
+		    tmap1 |= 1 << 0; // default set to the first bit: primary-vertex-used
+		}
+		else
+		{
+		    tmap1 = 0;
+		    tmap2 = 0;
+		    for (int iBit = 0; iBit < 32; iBit++)
+		    {
+			if (randy->Rndm() < 0.95)
+			tmap1 |= 1 << iBit;
+		    else 
+			tmap1 |= 0 << iBit;
+		    
+		    if(iBit < 27) 
+			if(randy->Rndm() < 0.95) 
+			    tmap2 |= 1 << iBit;
+			else 
+			    tmap2 |= 0 << iBit;
+		    }
+		    tmap1 |= 1 << 0; // default set to the first bit: primary-vertex-used
+		}
+	    }
+	    mHbtTrack->SetTopologyMap(0, tmap1);
+	    mHbtTrack->SetTopologyMap(1, tmap2);
+	
+	    mHbtTrack->SetTofBeta(ptot/p.Energy());
+	    mHbtTrack->SetHiddenInfo(0);
+	
+	    mHbtEvent->TrackCollection()->push_back(mHbtTrack);
+	} // end of a loop over tracks
+
+	//
+	// Number of charged tracks 
+	//
+	mHbtEvent->SetUncorrectedNumberOfPositivePrimaries(nTracksPos);
+	mHbtEvent->SetUncorrectedNumberOfNegativePrimaries(nTracksNeg);
+	
+    } // if event != 0
+    else
+    {
+	//
+	// FemtoDstMaker returned NULL, so we are not OK
+	//
+	mReaderStatus = 1;
+    }
+
+    return mHbtEvent;
+}
+
+//_________________
+StHbtEvent *StHbtStarGenHijingDstReader::ReturnHbtEvent() {
+
+    StHbtEvent *mHbtEvent = Read();
+    if (!mHbtEvent) {
+	std::cout << "StHbtStarGenHijingDstReader[WARN]: no hbt event (mHbtEvent == NULL)\n";
+	return 0;
+    }
+    else {
+	return mHbtEvent;
+    }
+}
+
+//_________________
+StHbtString StHbtStarGenHijingDstReader::Report() {
+
+    StHbtString temp = "\nThis is the StHbtStarGenHijingDstReader\n";
+    return temp;
+}
+
+//_________________
+void StHbtStarGenHijingDstReader::DoDebug(bool aDebug)
+{
+    mDebug = aDebug;
+}
+
+//________________
+float StHbtStarGenHijingDstReader::dedxMean(float mass, float momentum)
+{
+    Double_t dedxMean;
+    Double_t tpcDedxGain = 0.174325e-06;
+    Double_t tpcDedxOffset = -2.71889; 
+    Double_t tpcDedxRise = 776.626;
+  
+    Double_t gamma = TMath::Sqrt(momentum*momentum/(mass*mass)+1.);
+    Double_t beta = TMath::Sqrt(1. - 1./(gamma*gamma));
+    Double_t rise = tpcDedxRise* beta*beta*gamma*gamma;
+    
+    if (beta > 0)
+	dedxMean = tpcDedxGain/(beta*beta)*(0.5*TMath::Log(rise) - beta*beta 
+					    - tpcDedxOffset);
+    else
+	dedxMean = 1000.;
+    
+    return dedxMean;
+}

+ 77 - 0
StRoot/StFemtoDstMaker/StHbtStarGenHijingDstReader.h

@@ -0,0 +1,77 @@
+//
+//  The version of the class should be changed every time
+//  when any changes in the codes are done
+//  Grigory Nigmatkulov: 2016/12/15
+//
+#ifndef StHbtStarGenHijingDstReader_hh
+#define StHbtStarGenHijingDstReader_hh
+
+#include "StHbtMaker/Base/StHbtEventReader.hh"
+#include "StHbtMaker/Infrastructure/StHbtEvent.hh"
+#include "StHbtMaker/Infrastructure/StHbtEnumeration.hh"
+#include "StEvent/StEnumerations.h"
+
+#include "StMaker.h"
+#include "StFemtoEvent.h"
+#include "TSystem.h"
+#include "TChain.h"
+#include "StChain.h"
+#include "TRandom3.h"
+#include <StarGenerator/EVENT/StarGenParticle.h>
+#include <StarGenerator/EVENT/StarGenEvent.h>
+
+#include <string>
+#include <stdlib.h>
+
+class StFlowEvent;
+class StIOMaker;
+class TTree;
+class TFile;
+
+//_________________
+class StHbtStarGenHijingDstReader : public StHbtEventReader {
+private:
+    //
+    // General variables
+    //
+    string mDir;
+    string mFileName;
+    string mFilter;
+
+    //
+    // Internal variables
+    //
+    bool          mDebug;
+    unsigned int  mEventIndex;
+    unsigned int  mNEvents;
+    TChain       *mTChain;
+    TTree        *mTree;
+    StarGenEvent *mSGEvent;
+    float         mTotalTracks;
+    int           mMaxFiles;
+    TRandom3     *randy;
+
+    int InitRead(string dir, string fileName,
+		 string filter, int maxFiles);
+    int FillChain(TChain *aChain, const char *fileName, 
+		  int faxFiles);
+    int FillChain(TChain *chain, char *dir,
+		  const char *filter, int maxFiles);
+    void UninitRead();
+    float dedxMean(float mass, float momentum);
+    StHbtEvent *Read();
+  
+public:
+    StHbtStarGenHijingDstReader(const char *dirName, const char *fileName,
+				const char *filter = ".", int maxFiles = 1e9, bool debug = false);
+    StHbtEvent *ReturnHbtEvent();
+    void PrintTotalTracks();
+    void DoDebug(bool debug = false);
+    StHbtString Report();
+    int GetNEvents();
+  
+
+    ClassDef(StHbtStarGenHijingDstReader, 1)
+};
+
+#endif

+ 451 - 0
macros/HBT/expPionDCHBT_pau200_y2015_3D.C

@@ -0,0 +1,451 @@
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TString.h>
+#include <iostream>
+
+// #define MONITORS
+
+//Constants for analyses
+//Event information
+const int mPrimVertZbins = 6;                 //Number of z bins for analysis
+const double mPrimVertZ[2] = {-30., 30.};     //Redefine values of the vertex Z for each run [cm]
+const double mPrimVertR = 2.;                 //Distance in the XY plane [cm]
+const double mPrimVertVpdVz[2] = {-5., 5.};   //VpdZ - TPCZ difference
+const double mPrimVertXShift = 0.;
+const double mPrimVertYShift = 0.;
+const Bool_t mUseZdcCorrection = false;
+
+//Track information
+const double mChargedKaonMass = 0.493677;
+const double mChargedPionMass = 0.13957018;
+const int mPosCharge = +1;
+const int mNegCharge = -1;
+const int mTrackType = 1; 
+const int mPythiaKaonCode = 321;              //Pythia particle code
+const int mPythiaPionCode = 211;              //Pythia particle code
+
+const int mTrackNHits[2] = {15, 50};          //Number of track hits
+const double mTrackNHitsFit[2] = {15, 50};     //Number of track hits used in fit
+const double mTrackDCA[2] = {0., 2.};      //DCA of the track and PV [cm]
+const double mTrackASplitVal = 0.52;
+const double mTrackMom[2] = {0.15, 1.55};     //General momentum of the track [GeV/c]
+const double mTrackTransvMom[2] = {0.15, 1.55};  //General transverse momentum of the track [GeV/c]
+const double mTrackEta[2] = { -1., 1.};
+
+const short mDetectorSelection = 3;   
+
+const double mTrackTpcMom[2] = {0.15, 0.55};  //Track momentum for TPC idnetification [GeV/c]
+const double mnSigmaElectron[2] = {-1e9, 1e9};
+const double mnSigmaPion[2]     = {-2.,  2.};
+const double mnSigmaKaon[2]     = {-1e9, 1e9};
+const double mnSigmaProton[2]   = {-1e9, 1e9};
+
+const double mTrackTofMom[2] = {0.15, 1.55};   //Track momentum for ToF identification [GeV/c]
+const double mPionTofMsqr[2] = {-0.02, 0.062}; //Based on p+p data
+const double mKaonTofMsqr[2] = {0.20, 0.32};    //Based on p+p data
+
+const double mTrackTnTMom[2]   = {0.15, 1.55};
+const double mTnTNSigmaPion[2] = {-2., 2.};
+
+//Pair information
+const double mFmrCut[2] = {-1.1, 0.1};        //Not more than 10% of merged rows
+const double mSplitLevelCut[2] = {-0.5, 0.6}; //Splitting level from paper. StHbtPair::quality()
+const double mPairKt[2] = {0.,1.5};           //Pair transverse momentum cut
+const double mPairQinv[2] = {0., 1.5};        //Pair relative momentum
+const double mAveSep[2] = {5., 5000.};         //Average separation of tracks within TPC
+
+//Histograms values
+const int mPairKtBins = 10;
+const double mPairKtVal[2] = {0.1, 1.1};
+const int mQinvBins = 120;               // 10 MeV/c^2 bin width
+const double mQinvVal[2] = {0., 1.2};    // Redefine for each analysis
+
+const Char_t *defaultInFile = "sample.list";
+const Char_t *defaultOutFile = "oExpPion_pau200_y2015_3D.root";
+
+//_________________
+void expPionDCHBT_pau200_y2015_3D(const Char_t *inFileList = defaultInFile,
+				  const Char_t *oFileName = defaultOutFile,
+				  int mult = 0) {
+
+  std::cout << "****************************************" << std::endl
+	    << "*     expPionDCHBT_pau200_y2015_3D     *" << std::endl
+	    << "*                 Start                *" << std::endl
+	    << "****************************************" << std::endl << endl;
+
+  // load libraries
+  std::cout << "Loading libraries..." << std::endl;
+  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
+  loadSharedLibraries();
+  gSystem->Load("libMinuit");
+  gSystem->Load("StRefMultCorr");
+  gSystem->Load("StFlowMaker");   // should be placed before HBT
+  gSystem->Load("StHbtMaker");
+  gSystem->Load("StarClassLibrary");
+  gSystem->Load("StFemtoDstMaker");
+  gSystem->Load("libgsl");
+  gSystem->Load("libgslcblas");
+
+  std::cout << "Libraries have been successfully loaded" << std::endl;
+
+  // create chain
+  StChain* mChain = new StChain("StChain");
+  mChain->SetDebug(0);
+  StMuDebug::setLevel(0);
+
+  // set StHbtReader
+  std::cout << Form("Try to read file: %s", inFileList) << std::endl;
+  
+  StHbtFemtoDstReader *femtoReader = new StHbtFemtoDstReader("", inFileList, "root", -1, false);
+  /*
+                                    // trg             -- number of events
+  femtoReader->AddTrigId(500003);   // VPDMB-5-trgonly -- 24901523  
+  femtoReader->AddTrigId(500004);   // VPDMB-novtx     -- 23532384
+  femtoReader->AddTrigId(500904);   // VPDMB-30        -- 58381649
+                                    // total events    == 106815556
+  */
+  Long64_t mNEvents = femtoReader->GetNEvents();
+  
+  std::cout << Form("Number of events in TTree: %d", mNEvents) << std::endl;
+
+  // set StHbtMaker and StHbtManager
+  StHbtMaker *mHbtMaker = new StHbtMaker();
+  mHbtMaker->SetDebug(0);
+  StHbtManager *mHbtManager = mHbtMaker->HbtManager();
+  mHbtManager->SetEventReader(femtoReader);
+
+  // set basic event cut
+  gregsFemtoDstEventCut *mBasicEventCut = new gregsFemtoDstEventCut();
+  mBasicEventCut->SetCollisionType(0);     // 0-pp, 1-AuAu
+  mBasicEventCut->SetEventMult(0, 70);
+  mBasicEventCut->SetVertZPos(mPrimVertZ[0], mPrimVertZ[1]);
+  mBasicEventCut->SetVzVpdVzDiff(mPrimVertVpdVz[0],mPrimVertVpdVz[1]);
+  mBasicEventCut->SetVertRPos(mPrimVertR);
+  mBasicEventCut->SetEventCentralityBins(0, 8);  // for 5 centrality bins
+
+  // set basic pion single-particle cut
+  gregsPionTrackCut *mBasicTrackCut = new gregsPionTrackCut();
+  mBasicTrackCut->SetCharge(mPosCharge);
+  mBasicTrackCut->SetMass(mChargedPionMass);
+  mBasicTrackCut->SetType(mTrackType);
+  mBasicTrackCut->SetNHits(mTrackNHits[0], mTrackNHits[1]);
+  mBasicTrackCut->SetNHitsFit(mTrackNHitsFit[0], mTrackNHitsFit[1]);
+  mBasicTrackCut->SetAntiSplit(mTrackASplitVal);
+  mBasicTrackCut->SetEta(mTrackEta[0], mTrackEta[1]);
+  mBasicTrackCut->SetDCA(mTrackDCA[0], mTrackDCA[1]);
+  mBasicTrackCut->SetDCAGlobal(mTrackDCA[0], mTrackDCA[1]);
+  mBasicTrackCut->SetP(mTrackMom[0], mTrackMom[1]);
+
+  mBasicTrackCut->SetDetectorSelection(mDetectorSelection);
+
+  mBasicTrackCut->SetTpcMomentum(mTrackTpcMom[0], mTrackTpcMom[1]);
+  mBasicTrackCut->SetTpcNSigmaElectron(mnSigmaElectron[0],mnSigmaElectron[1]);
+  mBasicTrackCut->SetTpcNSigmaPion(mnSigmaPion[0], mnSigmaPion[1]);
+  mBasicTrackCut->SetTpcNSigmaKaon(mnSigmaKaon[0], mnSigmaKaon[1]);
+  mBasicTrackCut->SetTpcNSigmaProton(mnSigmaProton[0], mnSigmaProton[1]);
+
+  mBasicTrackCut->SetTofMassSqr(mPionTofMsqr[0],mPionTofMsqr[1]);
+  mBasicTrackCut->SetTofMomentum(mTrackTofMom[0], mTrackTofMom[1]);
+
+  mBasicTrackCut->SetTnTMomentum(mTrackTnTMom[0], mTrackTnTMom[1]);
+  mBasicTrackCut->SetTnTNSigmaPion(mTnTNSigmaPion[0], mTnTNSigmaPion[1]);
+
+  // set basic pion pair cut
+  gregsTrackPairCut *mBasicPairCut = new gregsTrackPairCut();
+  mBasicPairCut->SetFracOfMergedRow(mFmrCut[0],mFmrCut[1]);
+  mBasicPairCut->SetQuality(mSplitLevelCut[0],mSplitLevelCut[1]);
+  mBasicPairCut->SetKt(mPairKt[0],mPairKt[1]);
+  mBasicPairCut->SetQinv(mPairQinv[0],mPairQinv[1]);
+  mBasicPairCut->SetAverageSeparation(mAveSep[0],mAveSep[1]);
+
+  Int_t mLocalCharge, mMultBins, mMultLow, mMultHi, mEvents2Mix;
+  float ktLow, ktHi;
+
+  // The multiplicity analysis array, where
+  // the first element corresponds to the particle charge: (pos, neg),
+  // and 6 centralitiy bins: (1-100),(1-6),(7-9),(10-12),(13-16),(17,100)
+
+  // StHbtAnalysis, cut, CF and monitor declarations
+  StHbtVertexMultAnalysis *mPionMultAnalysis;
+  
+  // cuts
+  gregsFemtoDstEventCut *mPionMultEventCut;
+  gregsPionTrackCut     *mPionPosMultTrackCut, *mPionNegMultTrackCut;
+  gregsTrackPairCut     *mPionMultPairCut;
+
+  // correlation functions
+  QinvCorrFctnKt               *mPionQinvKtCF;
+  gregsBPLCMSFrame3DCorrFctnKt *mPionBPLCMSFrameKtCF;
+
+#ifdef MONITORS
+  //Cut monitors
+  fxtEventCutMonitor *mFxtPionEventCutMonitorPass;
+  fxtTrackCutMonitor *mFxtPionPosTrackCutMonitorPass, *mFxtPionNegTrackCutMonitorPass;
+  fxtPairCutMonitor  *mFxtPionPairCutMonitorPass;
+  fxtEventCutMonitor *mFxtPionEventCutMonitorFail;
+  fxtTrackCutMonitor *mFxtPionTrackCutMonitorFail;
+  fxtPairCutMonitor  *mFxtPionPairCutMonitorFail;
+#endif
+
+  std::cout << "---=== BEGIN CREATE ANALYSIS ===---\n";
+
+  switch(mult) {
+  case 1:  mMultLow = 0;  mMultHi = 6;    mMultBins = 2; mEvents2Mix = 5; break;
+  case 2:  mMultLow = 7;  mMultHi = 9;    mMultBins = 2; mEvents2Mix = 5; break;
+  case 3:  mMultLow = 10; mMultHi = 12;   mMultBins = 2; mEvents2Mix = 3; break;
+  case 4:  mMultLow = 13; mMultHi = 16;   mMultBins = 2; mEvents2Mix = 3; break;
+  case 5:  mMultLow = 17; mMultHi = 25;   mMultBins = 2; mEvents2Mix = 3; break;
+  case 6:  mMultLow = 26; mMultHi = 35;   mMultBins = 2; mEvents2Mix = 3; break;
+  case 7:  mMultLow = 36; mMultHi = 45;   mMultBins = 2; mEvents2Mix = 3; break;
+  case 8:  mMultLow = 46; mMultHi = 70;   mMultBins = 2; mEvents2Mix = 3; break;      
+  default: mMultLow = 0;  mMultHi = 70;   mMultBins = 7; mEvents2Mix = 3; break;
+  };
+
+  //Create analysis
+  mPionMultAnalysis = new StHbtVertexMultAnalysis(mPrimVertZbins,
+						  mPrimVertZ[0],
+						  mPrimVertZ[1],
+						  mMultBins,
+						  mMultLow,
+						  mMultHi);
+    
+  //Create event cut
+  mPionMultEventCut = new gregsFemtoDstEventCut(*mBasicEventCut);
+  mPionMultEventCut->SetEventMult(mMultLow, mMultHi);
+	
+#ifdef MONITORS
+  mFxtPionEventCutMonitorPass = 
+    new fxtEventCutMonitor("eventPionPass", Form("_%d", mult));
+  mFxtPionEventCutMonitorFail = 
+    new fxtEventCutMonitor("eventPionFail", Form("_%d", mult));
+  mPionEventCut->AddCutMonitor(mFxtPionEventCutMonitorPass,
+			       mFxtPionEventCutMonitorFail);
+#endif
+  mPionMultAnalysis->SetEventCut(mPionMultEventCut);
+
+  //
+  // Create single-track cuts and add them to the analysis
+  //
+  mPionPosMultTrackCut = new gregsPionTrackCut(*mBasicTrackCut);
+  mPionPosMultTrackCut->SetCharge(mPosCharge);
+  mPionNegMultTrackCut = new gregsPionTrackCut(*mBasicTrackCut);
+  mPionNegMultTrackCut->SetCharge(mNegCharge);
+    
+#ifdef MONITORS
+  //
+  // Positive charged pions cut monitors
+  //
+  mFxtPionPosTrackCutMonitorPass = 
+    new fxtTrackCutMonitor(Form("_trackPionPosPass_%d_", mult), mChargedPionMass);
+  mFxtPionPosTrackCutMonitorFail = 
+    new fxtTrackCutMonitor(Form("_trackPionPosFail_%d_", mult), mChargedPionMass);
+  mPionPosMultTrackCut->AddCutMonitor(mFxtPionPosTrackCutMonitorPass,
+				      mFxtPionPosTrackCutMonitorFail);
+  //
+  // Negative charged pions cut monitors
+  //
+  mFxtPionNegTrackCutMonitorPass = 
+    new fxtTrackCutMonitor(Form("_trackPionNegPass_%d_", mult), mChargedPionMass);
+  mFxtPionNegTrackCutMonitorFail = 
+    new fxtTrackCutMonitor(Form("_trackPionNegFail_%d_", mult), mChargedPionMass);
+  mPionNegMultTrackCut->AddCutMonitor(mFxtPionNegTrackCutMonitorPass,
+				      mFxtPionNegTrackCutMonitorFail);
+    
+#endif
+  mPionMultAnalysis->SetFirstParticleCut(mPionPosMultTrackCut);
+  mPionMultAnalysis->SetSecondParticleCut(mPionNegMultTrackCut);
+
+  //Create pair cut
+  mPionMultPairCut = new gregsTrackPairCut(*mBasicPairCut);
+
+#ifdef MONITORS
+  mFxtPionPairCutMonitorPass = 
+    new fxtPairCutMonitor(Form("_pairPionPass_%d_", mult));
+  mFxtPionPairCutMonitorPass->SetParticleMass(mChargedPionMass);
+  mFxtPionPairCutMonitorFail = 
+    new fxtPairCutMonitor(Form("_pairPionFail_%d_", mult));
+  mFxtPionPairCutMonitorFail->SetParticleMass(mChargedPionMass);
+  mPionPairCut->AddCutMonitor(mFxtPionPairCutMonitorPass,
+			      mFxtPionPairCutMonitorFail);
+#endif
+  mPionMultAnalysis->SetPairCut(mPionMultPairCut);
+  //
+  // Correlation function initialization
+  //
+  // 1D
+  mPionQinvKtCF = 
+    new QinvCorrFctnKt(Form("hPionDC_QinvKtMix_%d", mult),
+		       mQinvBins, mQinvVal[0], mQinvVal[1],
+		       mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+  // 3D
+  mPionBPLCMSFrameKtCF = 
+    new gregsBPLCMSFrame3DCorrFctnKt(Form("hPionDC_BPLCMSFrameKtCF_%d", mult),
+				     mQinvBins, mQinvVal[0], mQinvVal[1],
+				     mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+
+  //
+  // Add correlation functions to the HbtManager
+  //
+  mPionMultAnalysis->AddCorrFctn(mPionQinvKtCF); 
+  mPionMultAnalysis->AddCorrFctn(mPionBPLCMSFrameKtCF);   
+  mPionMultAnalysis->SetNumEventsToMix(mEvents2Mix);   // set number of events to mix
+  mHbtManager->AddAnalysis(mPionMultAnalysis); // add analysis to the manager
+  
+  std::cout << "---=== END CREATE ANALYSIS ===---\n";
+  
+  mChain->Init(); // StChain initialization
+
+  //
+  // Main loop over all events
+  //
+  Int_t iRet = 0;
+  Int_t mNEventsProcessed = 0;
+  Float_t mPercentCounter = 0.05;
+  Float_t mProgress = 0.;
+
+  unsigned int mSubCounter = 1;              // time
+  time_t mStartTime, mStopTime, mDiffTime;   // time
+  float mFrac;                               // time
+  mStartTime = time(0);                      // time
+
+  for(Int_t iEvent=0; iEvent<mNEvents; iEvent++) {
+
+    mProgress = (Float_t)iEvent / (Float_t)mNEvents;
+    mNEventsProcessed++;
+    
+    if(mProgress >= mPercentCounter) {
+      mPercentCounter += 0.05;
+      mStopTime = time(0);
+      mDiffTime = difftime(mStopTime, mStartTime);
+      mFrac = (float)mDiffTime*(float)(mNEvents - iEvent)/(float)iEvent;
+      mSubCounter++;
+      std::cout << Form("Processing progress: %4.2f%%. Time left: %.1f sec", mProgress*100., mFrac)
+		<< std::endl;
+    }
+
+    mChain->Clear();
+    iRet = mChain->Make(iEvent);
+    if(iRet!=0) {
+      std::cout << "expPionHBT_pau200_y2015_3D[ERROR] - iRet: " << iRet << std::endl;
+      break;
+    }
+  }
+
+  mChain->Finish();
+  
+  TFile *oFile = new TFile(oFileName, "RECREATE"); // create output file
+  std::cout << "Writing histograms to the output file...";
+
+  
+#ifdef MONITORS
+  //
+  // Event cut monitors
+  //
+  mFxtPionEventCutMonitorPass->VertexZ()->Write();
+  mFxtPionEventCutMonitorPass->VertexYvsVertexX()->Write();
+  mFxtPionEventCutMonitorPass->RefMult()->Write();
+  mFxtPionEventCutMonitorPass->VpdVzDiff()->Write();
+  mFxtPionEventCutMonitorFail->VertexZ()->Write();
+  mFxtPionEventCutMonitorFail->VertexYvsVertexX()->Write();
+  mFxtPionEventCutMonitorFail->RefMult()->Write();
+  mFxtPionEventCutMonitorFail->VpdVzDiff()->Write();
+  //
+  // Positive pion track cut monitors
+  //
+  mFxtPionPosTrackCutMonitorPass->DCAGlobal()->Write();
+  mFxtPionPosTrackCutMonitorPass->Nhits()->Write();
+  mFxtPionPosTrackCutMonitorPass->P()->Write();
+  mFxtPionPosTrackCutMonitorPass->Pt()->Write();
+  mFxtPionPosTrackCutMonitorPass->PtVsNsigmaPion()->Write();
+  mFxtPionPosTrackCutMonitorPass->PtVsNsigmaKaon()->Write();
+  mFxtPionPosTrackCutMonitorPass->PtVsNsigmaProton()->Write();
+  mFxtPionPosTrackCutMonitorPass->PvsDedx()->Write();
+  mFxtPionPosTrackCutMonitorPass->Rapidity()->Write();
+  mFxtPionPosTrackCutMonitorPass->PseudoRapidity()->Write();
+  mFxtPionPosTrackCutMonitorPass->PtVsMassSqr()->Write();
+  mFxtPionPosTrackCutMonitorFail->DCAGlobal()->Write();
+  mFxtPionPosTrackCutMonitorFail->Nhits()->Write();
+  mFxtPionPosTrackCutMonitorFail->P()->Write();
+  mFxtPionPosTrackCutMonitorFail->Pt()->Write();
+  mFxtPionPosTrackCutMonitorFail->PtVsNsigmaPion()->Write();
+  mFxtPionPosTrackCutMonitorFail->PtVsNsigmaKaon()->Write();
+  mFxtPionPosTrackCutMonitorFail->PtVsNsigmaProton()->Write();
+  mFxtPionPosTrackCutMonitorFail->PvsDedx()->Write();
+  mFxtPionPosTrackCutMonitorFail->Rapidity()->Write();
+  mFxtPionPosTrackCutMonitorFail->PseudoRapidity()->Write();
+  mFxtPionPosTrackCutMonitorFail->PtVsMassSqr()->Write();
+  //
+  // Negative pion track cut monitors
+  //
+  mFxtPionNegTrackCutMonitorPass->DCAGlobal()->Write();
+  mFxtPionNegTrackCutMonitorPass->Nhits()->Write();
+  mFxtPionNegTrackCutMonitorPass->P()->Write();
+  mFxtPionNegTrackCutMonitorPass->Pt()->Write();
+  mFxtPionNegTrackCutMonitorPass->PtVsNsigmaPion()->Write();
+  mFxtPionNegTrackCutMonitorPass->PtVsNsigmaKaon()->Write();
+  mFxtPionNegTrackCutMonitorPass->PtVsNsigmaProton()->Write();
+  mFxtPionNegTrackCutMonitorPass->PvsDedx()->Write();
+  mFxtPionNegTrackCutMonitorPass->Rapidity()->Write();
+  mFxtPionNegTrackCutMonitorPass->PseudoRapidity()->Write();
+  mFxtPionNegTrackCutMonitorPass->PtVsMassSqr()->Write();
+  mFxtPionNegTrackCutMonitorFail->DCAGlobal()->Write();
+  mFxtPionNegTrackCutMonitorFail->Nhits()->Write();
+  mFxtPionNegTrackCutMonitorFail->P()->Write();
+  mFxtPionNegTrackCutMonitorFail->Pt()->Write();
+  mFxtPionNegTrackCutMonitorFail->PtVsNsigmaPion()->Write();
+  mFxtPionNegTrackCutMonitorFail->PtVsNsigmaKaon()->Write();
+  mFxtPionNegTrackCutMonitorFail->PtVsNsigmaProton()->Write();
+  mFxtPionNegTrackCutMonitorFail->PvsDedx()->Write();
+  mFxtPionNegTrackCutMonitorFail->Rapidity()->Write();
+  mFxtPionNegTrackCutMonitorFail->PseudoRapidity()->Write();
+  mFxtPionNegTrackCutMonitorFail->PtVsMassSqr()->Write();
+  //
+  // Pion pair cut monitors
+  //
+  mFxtPionPairCutMonitorPass->Kt()->Write();
+  mFxtPionPairCutMonitorPass->Pt1Pt2Qinv()->Write();
+  mFxtPionPairCutMonitorPass->FMRvsQinv()->Write();
+  mFxtPionPairCutMonitorPass->SLvsQinv()->Write();
+  mFxtPionPairCutMonitorPass->AveSepVsQinv()->Write();
+  mFxtPionPairCutMonitorFail->Kt()->Write();
+  mFxtPionPairCutMonitorFail->Pt1Pt2Qinv()->Write();
+  mFxtPionPairCutMonitorFail->FMRvsQinv()->Write();
+  mFxtPionPairCutMonitorFail->SLvsQinv()->Write();
+  mFxtPionPairCutMonitorFail->AveSepVsQinv()->Write();
+#endif
+
+  //
+  // Write correlation function to the
+  //   output file
+  //
+  for (int iKt = 0; iKt < mPairKtBins; iKt++) {
+    mPionQinvKtCF->Numerator(iKt)->Write();
+    mPionQinvKtCF->Denominator(iKt)->Write();
+      
+    mPionBPLCMSFrameKtCF->Numerator(iKt)->Write();
+    mPionBPLCMSFrameKtCF->Denominator(iKt)->Write();
+    mPionBPLCMSFrameKtCF->QinvHisto(iKt)->Write();
+  }
+  std::cout << "\t[DONE]" << std::endl;
+
+  //
+  // Write the output file and close it
+  //
+  std::cout << "Saving the output file...";
+  oFile->Write();
+  oFile->Close();
+  std::cout << "\t[DONE]" << std::endl;
+
+  //
+  // Delete chain and prompt
+  //   last message :D
+  //
+  delete mChain;
+
+  std::cout << "****************************************" << std::endl
+	    << "*      expPionDCHBT_pau200_y2015_3D    *" << std::endl
+	    << "*                Finish                *" << std::endl
+	    << "****************************************" << std::endl << std::endl;
+}

+ 444 - 0
macros/HBT/simPionDCHBT_pau200_3D.C

@@ -0,0 +1,444 @@
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TString.h>
+#include <iostream>
+
+// #define MONITORS
+
+//Constants for analyses
+//Event information
+const int mPrimVertZbins = 6;                 //Number of z bins for analysis
+const double mPrimVertZ[2] = {-30., 30.};     //Redefine values of the vertex Z for each run [cm]
+const double mPrimVertR = 2.;                 //Distance in the XY plane [cm]
+const double mPrimVertVpdVz[2] = {-5., 5.};   //VpdZ - TPCZ difference
+const double mPrimVertXShift = 0.;
+const double mPrimVertYShift = 0.;
+const Bool_t mUseZdcCorrection = false;
+
+//Track information
+const double mChargedKaonMass = 0.493677;
+const double mChargedPionMass = 0.13957018;
+const int mPosCharge = +1;
+const int mNegCharge = -1;
+const int mTrackType = 1; 
+const int mPythiaKaonCode = 321;              //Pythia particle code
+const int mPythiaPionCode = 211;              //Pythia particle code
+
+const int mTrackNHits[2] = {15, 50};          //Number of track hits
+const double mTrackNHitsFit[2] = {15, 50};     //Number of track hits used in fit
+const double mTrackDCA[2] = {0., 2.};      //DCA of the track and PV [cm]
+const double mTrackASplitVal = 0.52;
+const double mTrackMom[2] = {0.15, 1.55};     //General momentum of the track [GeV/c]
+const double mTrackTransvMom[2] = {0.15, 1.55};  //General transverse momentum of the track [GeV/c]
+const double mTrackEta[2] = { -1., 1.};
+
+const short mDetectorSelection = 3;   
+
+const double mTrackTpcMom[2] = {0.15, 0.55};  //Track momentum for TPC idnetification [GeV/c]
+const double mnSigmaElectron[2] = {-1e9, 1e9};
+const double mnSigmaPion[2]     = {-2.,  2.};
+const double mnSigmaKaon[2]     = {-1e9, 1e9};
+const double mnSigmaProton[2]   = {-1e9, 1e9};
+
+const double mTrackTofMom[2] = {0.15, 1.55};   //Track momentum for ToF identification [GeV/c]
+const double mPionTofMsqr[2] = {-0.02, 0.062}; //Based on p+p data
+const double mKaonTofMsqr[2] = {0.20, 0.32};    //Based on p+p data
+
+const double mTrackTnTMom[2]   = {0.15, 1.55};
+const double mTnTNSigmaPion[2] = {-2., 2.};
+
+//Pair information
+const double mFmrCut[2] = {-1.1, 0.1};        //Not more than 10% of merged rows
+const double mSplitLevelCut[2] = {-0.5, 0.6}; //Splitting level from paper. StHbtPair::quality()
+const double mPairKt[2] = {0.,1.5};           //Pair transverse momentum cut
+const double mPairQinv[2] = {0., 1.5};        //Pair relative momentum
+const double mAveSep[2] = {5., 5000.};         //Average separation of tracks within TPC
+
+//Histograms values
+const int mPairKtBins = 10;
+const double mPairKtVal[2] = {0.1, 1.1};
+const int mQinvBins = 120;               // 10 MeV/c^2 bin width
+const double mQinvVal[2] = {0., 1.2};    // Redefine for each analysis
+
+const Char_t *defaultInFile = "sample.list";
+const Char_t *defaultOutFile = "oSimPion_pau200_3D.root";
+
+//_________________
+void expPionDCHBT_pau200_y2015_3D(const Char_t *inFileList = defaultInFile,
+				  const Char_t *oFileName = defaultOutFile,
+				  int mult = 0) {
+
+  std::cout << "****************************************" << std::endl
+	    << "*     simPionDCHBT_pau200_3D     *" << std::endl
+	    << "*                 Start                *" << std::endl
+	    << "****************************************" << std::endl << endl;
+
+  // load libraries
+  std::cout << "Loading libraries..." << std::endl;
+  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
+  loadSharedLibraries();
+  gSystem->Load("libMinuit");
+  gSystem->Load("StRefMultCorr");
+  gSystem->Load("StFlowMaker");   // should be placed before HBT
+  gSystem->Load("StHbtMaker");
+  gSystem->Load("StarClassLibrary");
+  gSystem->Load("StFemtoDstMaker");
+  gSystem->Load("libgsl");
+  gSystem->Load("libgslcblas");
+
+  std::cout << "Libraries have been successfully loaded" << std::endl;
+
+  // create chain
+  StChain* mChain = new StChain("StChain");
+  mChain->SetDebug(0);
+  StMuDebug::setLevel(0);
+
+  // set StHbtReader
+  std::cout << Form("Try to read file: %s", inFileList) << std::endl;
+  
+  StHbtStarGenHijingDstReader *hijingReader = new StHbtStarGenHijingDstReader("", inFileList, "root", -1, false);
+  Long64_t mNEvents = hijingReader->GetNEvents();
+  
+  std::cout << Form("Number of events in TTree: %d", mNEvents) << std::endl;
+
+  // set StHbtMaker and StHbtManager
+  StHbtMaker *mHbtMaker = new StHbtMaker();
+  mHbtMaker->SetDebug(0);
+  StHbtManager *mHbtManager = mHbtMaker->HbtManager();
+  mHbtManager->SetEventReader(hijingReader);
+
+  // set basic event cut
+  gregsFemtoDstEventCut *mBasicEventCut = new gregsFemtoDstEventCut();
+  mBasicEventCut->SetCollisionType(0);     // 0-pp, 1-AuAu
+  mBasicEventCut->SetEventMult(0, 70);
+  mBasicEventCut->SetVertZPos(mPrimVertZ[0], mPrimVertZ[1]);
+  mBasicEventCut->SetVzVpdVzDiff(mPrimVertVpdVz[0],mPrimVertVpdVz[1]);
+  mBasicEventCut->SetVertRPos(mPrimVertR);
+  mBasicEventCut->SetEventCentralityBins(0, 8);  // for 5 centrality bins
+
+  // set basic pion single-particle cut
+  gregsPionTrackCut *mBasicTrackCut = new gregsPionTrackCut();
+  mBasicTrackCut->SetCharge(mPosCharge);
+  mBasicTrackCut->SetMass(mChargedPionMass);
+  mBasicTrackCut->SetType(mTrackType);
+  mBasicTrackCut->SetNHits(mTrackNHits[0], mTrackNHits[1]);
+  mBasicTrackCut->SetNHitsFit(mTrackNHitsFit[0], mTrackNHitsFit[1]);
+  mBasicTrackCut->SetAntiSplit(mTrackASplitVal);
+  mBasicTrackCut->SetEta(mTrackEta[0], mTrackEta[1]);
+  mBasicTrackCut->SetDCA(mTrackDCA[0], mTrackDCA[1]);
+  mBasicTrackCut->SetDCAGlobal(mTrackDCA[0], mTrackDCA[1]);
+  mBasicTrackCut->SetP(mTrackMom[0], mTrackMom[1]);
+
+  mBasicTrackCut->SetDetectorSelection(mDetectorSelection);
+
+  mBasicTrackCut->SetTpcMomentum(mTrackTpcMom[0], mTrackTpcMom[1]);
+  mBasicTrackCut->SetTpcNSigmaElectron(mnSigmaElectron[0],mnSigmaElectron[1]);
+  mBasicTrackCut->SetTpcNSigmaPion(mnSigmaPion[0], mnSigmaPion[1]);
+  mBasicTrackCut->SetTpcNSigmaKaon(mnSigmaKaon[0], mnSigmaKaon[1]);
+  mBasicTrackCut->SetTpcNSigmaProton(mnSigmaProton[0], mnSigmaProton[1]);
+
+  mBasicTrackCut->SetTofMassSqr(mPionTofMsqr[0],mPionTofMsqr[1]);
+  mBasicTrackCut->SetTofMomentum(mTrackTofMom[0], mTrackTofMom[1]);
+
+  mBasicTrackCut->SetTnTMomentum(mTrackTnTMom[0], mTrackTnTMom[1]);
+  mBasicTrackCut->SetTnTNSigmaPion(mTnTNSigmaPion[0], mTnTNSigmaPion[1]);
+
+  // set basic pion pair cut
+  gregsTrackPairCut *mBasicPairCut = new gregsTrackPairCut();
+  mBasicPairCut->SetFracOfMergedRow(mFmrCut[0],mFmrCut[1]);
+  mBasicPairCut->SetQuality(mSplitLevelCut[0],mSplitLevelCut[1]);
+  mBasicPairCut->SetKt(mPairKt[0],mPairKt[1]);
+  mBasicPairCut->SetQinv(mPairQinv[0],mPairQinv[1]);
+  mBasicPairCut->SetAverageSeparation(mAveSep[0],mAveSep[1]);
+
+  Int_t mLocalCharge, mMultBins, mMultLow, mMultHi, mEvents2Mix;
+  float ktLow, ktHi;
+
+  // The multiplicity analysis array, where
+  // the first element corresponds to the particle charge: (pos, neg),
+  // and 6 centralitiy bins: (1-100),(1-6),(7-9),(10-12),(13-16),(17,100)
+
+  // StHbtAnalysis, cut, CF and monitor declarations
+  StHbtVertexMultAnalysis *mPionMultAnalysis;
+  
+  // cuts
+  gregsFemtoDstEventCut *mPionMultEventCut;
+  gregsPionTrackCut     *mPionPosMultTrackCut, *mPionNegMultTrackCut;
+  gregsTrackPairCut     *mPionMultPairCut;
+
+  // correlation functions
+  QinvCorrFctnKt               *mPionQinvKtCF;
+  gregsBPLCMSFrame3DCorrFctnKt *mPionBPLCMSFrameKtCF;
+
+#ifdef MONITORS
+  //Cut monitors
+  fxtEventCutMonitor *mFxtPionEventCutMonitorPass;
+  fxtTrackCutMonitor *mFxtPionPosTrackCutMonitorPass, *mFxtPionNegTrackCutMonitorPass;
+  fxtPairCutMonitor  *mFxtPionPairCutMonitorPass;
+  fxtEventCutMonitor *mFxtPionEventCutMonitorFail;
+  fxtTrackCutMonitor *mFxtPionTrackCutMonitorFail;
+  fxtPairCutMonitor  *mFxtPionPairCutMonitorFail;
+#endif
+
+  std::cout << "---=== BEGIN CREATE ANALYSIS ===---\n";
+
+  switch(mult) {
+  case 1:  mMultLow = 0;  mMultHi = 6;    mMultBins = 2; mEvents2Mix = 5; break;
+  case 2:  mMultLow = 7;  mMultHi = 9;    mMultBins = 2; mEvents2Mix = 5; break;
+  case 3:  mMultLow = 10; mMultHi = 12;   mMultBins = 2; mEvents2Mix = 3; break;
+  case 4:  mMultLow = 13; mMultHi = 16;   mMultBins = 2; mEvents2Mix = 3; break;
+  case 5:  mMultLow = 17; mMultHi = 25;   mMultBins = 2; mEvents2Mix = 3; break;
+  case 6:  mMultLow = 26; mMultHi = 35;   mMultBins = 2; mEvents2Mix = 3; break;
+  case 7:  mMultLow = 36; mMultHi = 45;   mMultBins = 2; mEvents2Mix = 3; break;
+  case 8:  mMultLow = 46; mMultHi = 70;   mMultBins = 2; mEvents2Mix = 3; break;      
+  default: mMultLow = 0;  mMultHi = 70;   mMultBins = 7; mEvents2Mix = 3; break;
+  };
+
+  //Create analysis
+  mPionMultAnalysis = new StHbtVertexMultAnalysis(mPrimVertZbins,
+						  mPrimVertZ[0],
+						  mPrimVertZ[1],
+						  mMultBins,
+						  mMultLow,
+						  mMultHi);
+    
+  //Create event cut
+  mPionMultEventCut = new gregsFemtoDstEventCut(*mBasicEventCut);
+  mPionMultEventCut->SetEventMult(mMultLow, mMultHi);
+	
+#ifdef MONITORS
+  mFxtPionEventCutMonitorPass = 
+    new fxtEventCutMonitor("eventPionPass", Form("_%d", mult));
+  mFxtPionEventCutMonitorFail = 
+    new fxtEventCutMonitor("eventPionFail", Form("_%d", mult));
+  mPionEventCut->AddCutMonitor(mFxtPionEventCutMonitorPass,
+			       mFxtPionEventCutMonitorFail);
+#endif
+  mPionMultAnalysis->SetEventCut(mPionMultEventCut);
+
+  //
+  // Create single-track cuts and add them to the analysis
+  //
+  mPionPosMultTrackCut = new gregsPionTrackCut(*mBasicTrackCut);
+  mPionPosMultTrackCut->SetCharge(mPosCharge);
+  mPionNegMultTrackCut = new gregsPionTrackCut(*mBasicTrackCut);
+  mPionNegMultTrackCut->SetCharge(mNegCharge);
+    
+#ifdef MONITORS
+  //
+  // Positive charged pions cut monitors
+  //
+  mFxtPionPosTrackCutMonitorPass = 
+    new fxtTrackCutMonitor(Form("_trackPionPosPass_%d_", mult), mChargedPionMass);
+  mFxtPionPosTrackCutMonitorFail = 
+    new fxtTrackCutMonitor(Form("_trackPionPosFail_%d_", mult), mChargedPionMass);
+  mPionPosMultTrackCut->AddCutMonitor(mFxtPionPosTrackCutMonitorPass,
+				      mFxtPionPosTrackCutMonitorFail);
+  //
+  // Negative charged pions cut monitors
+  //
+  mFxtPionNegTrackCutMonitorPass = 
+    new fxtTrackCutMonitor(Form("_trackPionNegPass_%d_", mult), mChargedPionMass);
+  mFxtPionNegTrackCutMonitorFail = 
+    new fxtTrackCutMonitor(Form("_trackPionNegFail_%d_", mult), mChargedPionMass);
+  mPionNegMultTrackCut->AddCutMonitor(mFxtPionNegTrackCutMonitorPass,
+				      mFxtPionNegTrackCutMonitorFail);
+    
+#endif
+  mPionMultAnalysis->SetFirstParticleCut(mPionPosMultTrackCut);
+  mPionMultAnalysis->SetSecondParticleCut(mPionNegMultTrackCut);
+
+  //Create pair cut
+  mPionMultPairCut = new gregsTrackPairCut(*mBasicPairCut);
+
+#ifdef MONITORS
+  mFxtPionPairCutMonitorPass = 
+    new fxtPairCutMonitor(Form("_pairPionPass_%d_", mult));
+  mFxtPionPairCutMonitorPass->SetParticleMass(mChargedPionMass);
+  mFxtPionPairCutMonitorFail = 
+    new fxtPairCutMonitor(Form("_pairPionFail_%d_", mult));
+  mFxtPionPairCutMonitorFail->SetParticleMass(mChargedPionMass);
+  mPionPairCut->AddCutMonitor(mFxtPionPairCutMonitorPass,
+			      mFxtPionPairCutMonitorFail);
+#endif
+  mPionMultAnalysis->SetPairCut(mPionMultPairCut);
+  //
+  // Correlation function initialization
+  //
+  // 1D
+  mPionQinvKtCF = 
+    new QinvCorrFctnKt(Form("hPionDC_QinvKtMix_%d", mult),
+		       mQinvBins, mQinvVal[0], mQinvVal[1],
+		       mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+  // 3D
+  mPionBPLCMSFrameKtCF = 
+    new gregsBPLCMSFrame3DCorrFctnKt(Form("hPionDC_BPLCMSFrameKtCF_%d", mult),
+				     mQinvBins, mQinvVal[0], mQinvVal[1],
+				     mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+
+  //
+  // Add correlation functions to the HbtManager
+  //
+  mPionMultAnalysis->AddCorrFctn(mPionQinvKtCF); 
+  mPionMultAnalysis->AddCorrFctn(mPionBPLCMSFrameKtCF);   
+  mPionMultAnalysis->SetNumEventsToMix(mEvents2Mix);   // set number of events to mix
+  mHbtManager->AddAnalysis(mPionMultAnalysis); // add analysis to the manager
+  
+  std::cout << "---=== END CREATE ANALYSIS ===---\n";
+  
+  mChain->Init(); // StChain initialization
+
+  //
+  // Main loop over all events
+  //
+  Int_t iRet = 0;
+  Int_t mNEventsProcessed = 0;
+  Float_t mPercentCounter = 0.05;
+  Float_t mProgress = 0.;
+
+  unsigned int mSubCounter = 1;              // time
+  time_t mStartTime, mStopTime, mDiffTime;   // time
+  float mFrac;                               // time
+  mStartTime = time(0);                      // time
+
+  for(Int_t iEvent=0; iEvent<mNEvents; iEvent++) {
+
+    mProgress = (Float_t)iEvent / (Float_t)mNEvents;
+    mNEventsProcessed++;
+    
+    if(mProgress >= mPercentCounter) {
+      mPercentCounter += 0.05;
+      mStopTime = time(0);
+      mDiffTime = difftime(mStopTime, mStartTime);
+      mFrac = (float)mDiffTime*(float)(mNEvents - iEvent)/(float)iEvent;
+      mSubCounter++;
+      std::cout << Form("Processing progress: %4.2f%%. Time left: %.1f sec", mProgress*100., mFrac)
+		<< std::endl;
+    }
+
+    mChain->Clear();
+    iRet = mChain->Make(iEvent);
+    if(iRet!=0) {
+      std::cout << "expPionHBT_pau200_y2015_3D[ERROR] - iRet: " << iRet << std::endl;
+      break;
+    }
+  }
+
+  mChain->Finish();
+  
+  TFile *oFile = new TFile(oFileName, "RECREATE"); // create output file
+  std::cout << "Writing histograms to the output file...";
+
+  
+#ifdef MONITORS
+  //
+  // Event cut monitors
+  //
+  mFxtPionEventCutMonitorPass->VertexZ()->Write();
+  mFxtPionEventCutMonitorPass->VertexYvsVertexX()->Write();
+  mFxtPionEventCutMonitorPass->RefMult()->Write();
+  mFxtPionEventCutMonitorPass->VpdVzDiff()->Write();
+  mFxtPionEventCutMonitorFail->VertexZ()->Write();
+  mFxtPionEventCutMonitorFail->VertexYvsVertexX()->Write();
+  mFxtPionEventCutMonitorFail->RefMult()->Write();
+  mFxtPionEventCutMonitorFail->VpdVzDiff()->Write();
+  //
+  // Positive pion track cut monitors
+  //
+  mFxtPionPosTrackCutMonitorPass->DCAGlobal()->Write();
+  mFxtPionPosTrackCutMonitorPass->Nhits()->Write();
+  mFxtPionPosTrackCutMonitorPass->P()->Write();
+  mFxtPionPosTrackCutMonitorPass->Pt()->Write();
+  mFxtPionPosTrackCutMonitorPass->PtVsNsigmaPion()->Write();
+  mFxtPionPosTrackCutMonitorPass->PtVsNsigmaKaon()->Write();
+  mFxtPionPosTrackCutMonitorPass->PtVsNsigmaProton()->Write();
+  mFxtPionPosTrackCutMonitorPass->PvsDedx()->Write();
+  mFxtPionPosTrackCutMonitorPass->Rapidity()->Write();
+  mFxtPionPosTrackCutMonitorPass->PseudoRapidity()->Write();
+  mFxtPionPosTrackCutMonitorPass->PtVsMassSqr()->Write();
+  mFxtPionPosTrackCutMonitorFail->DCAGlobal()->Write();
+  mFxtPionPosTrackCutMonitorFail->Nhits()->Write();
+  mFxtPionPosTrackCutMonitorFail->P()->Write();
+  mFxtPionPosTrackCutMonitorFail->Pt()->Write();
+  mFxtPionPosTrackCutMonitorFail->PtVsNsigmaPion()->Write();
+  mFxtPionPosTrackCutMonitorFail->PtVsNsigmaKaon()->Write();
+  mFxtPionPosTrackCutMonitorFail->PtVsNsigmaProton()->Write();
+  mFxtPionPosTrackCutMonitorFail->PvsDedx()->Write();
+  mFxtPionPosTrackCutMonitorFail->Rapidity()->Write();
+  mFxtPionPosTrackCutMonitorFail->PseudoRapidity()->Write();
+  mFxtPionPosTrackCutMonitorFail->PtVsMassSqr()->Write();
+  //
+  // Negative pion track cut monitors
+  //
+  mFxtPionNegTrackCutMonitorPass->DCAGlobal()->Write();
+  mFxtPionNegTrackCutMonitorPass->Nhits()->Write();
+  mFxtPionNegTrackCutMonitorPass->P()->Write();
+  mFxtPionNegTrackCutMonitorPass->Pt()->Write();
+  mFxtPionNegTrackCutMonitorPass->PtVsNsigmaPion()->Write();
+  mFxtPionNegTrackCutMonitorPass->PtVsNsigmaKaon()->Write();
+  mFxtPionNegTrackCutMonitorPass->PtVsNsigmaProton()->Write();
+  mFxtPionNegTrackCutMonitorPass->PvsDedx()->Write();
+  mFxtPionNegTrackCutMonitorPass->Rapidity()->Write();
+  mFxtPionNegTrackCutMonitorPass->PseudoRapidity()->Write();
+  mFxtPionNegTrackCutMonitorPass->PtVsMassSqr()->Write();
+  mFxtPionNegTrackCutMonitorFail->DCAGlobal()->Write();
+  mFxtPionNegTrackCutMonitorFail->Nhits()->Write();
+  mFxtPionNegTrackCutMonitorFail->P()->Write();
+  mFxtPionNegTrackCutMonitorFail->Pt()->Write();
+  mFxtPionNegTrackCutMonitorFail->PtVsNsigmaPion()->Write();
+  mFxtPionNegTrackCutMonitorFail->PtVsNsigmaKaon()->Write();
+  mFxtPionNegTrackCutMonitorFail->PtVsNsigmaProton()->Write();
+  mFxtPionNegTrackCutMonitorFail->PvsDedx()->Write();
+  mFxtPionNegTrackCutMonitorFail->Rapidity()->Write();
+  mFxtPionNegTrackCutMonitorFail->PseudoRapidity()->Write();
+  mFxtPionNegTrackCutMonitorFail->PtVsMassSqr()->Write();
+  //
+  // Pion pair cut monitors
+  //
+  mFxtPionPairCutMonitorPass->Kt()->Write();
+  mFxtPionPairCutMonitorPass->Pt1Pt2Qinv()->Write();
+  mFxtPionPairCutMonitorPass->FMRvsQinv()->Write();
+  mFxtPionPairCutMonitorPass->SLvsQinv()->Write();
+  mFxtPionPairCutMonitorPass->AveSepVsQinv()->Write();
+  mFxtPionPairCutMonitorFail->Kt()->Write();
+  mFxtPionPairCutMonitorFail->Pt1Pt2Qinv()->Write();
+  mFxtPionPairCutMonitorFail->FMRvsQinv()->Write();
+  mFxtPionPairCutMonitorFail->SLvsQinv()->Write();
+  mFxtPionPairCutMonitorFail->AveSepVsQinv()->Write();
+#endif
+
+  //
+  // Write correlation function to the
+  //   output file
+  //
+  for (int iKt = 0; iKt < mPairKtBins; iKt++) {
+    mPionQinvKtCF->Numerator(iKt)->Write();
+    mPionQinvKtCF->Denominator(iKt)->Write();
+      
+    mPionBPLCMSFrameKtCF->Numerator(iKt)->Write();
+    mPionBPLCMSFrameKtCF->Denominator(iKt)->Write();
+    mPionBPLCMSFrameKtCF->QinvHisto(iKt)->Write();
+  }
+  std::cout << "\t[DONE]" << std::endl;
+
+  //
+  // Write the output file and close it
+  //
+  std::cout << "Saving the output file...";
+  oFile->Write();
+  oFile->Close();
+  std::cout << "\t[DONE]" << std::endl;
+
+  //
+  // Delete chain and prompt
+  //   last message :D
+  //
+  delete mChain;
+
+  std::cout << "****************************************" << std::endl
+	    << "*      expPionDCHBT_pau200_y2015_3D    *" << std::endl
+	    << "*                Finish                *" << std::endl
+	    << "****************************************" << std::endl << std::endl;
+}

+ 194 - 0
macros/femtoDst/FemtoDstMaker_mesons_pp500_y2011.C

@@ -0,0 +1,194 @@
+#include "TROOT.h"
+#include "TSystem.h"
+#include "TChain.h"
+#include "TFile.h"
+#include "TString.h"
+#include <iostream>
+
+const Char_t *defaultInFile  = "/star/data01/pwg/pusheax/data/test/pp500/st_physics_12083058_raw_3020003.MuDst.root";
+const Char_t *defaultOutFile = "oFemtoDstMesons_pp500_2011.root";
+
+class StFemtoDstInclusiveSelector;
+StFemtoDstInclusiveSelector *mFemtoDstInclSelector = NULL;
+
+//_________________
+void FemtoDstMaker_mesons_pp500_y2011(const Char_t* inFileList = defaultInFile,
+				      const Char_t* outFileName = defaultOutFile,
+				      const Bool_t  isAuAu = false) {
+  
+  std::cout << "**************************" << std::endl
+	    << "*       FemtoDstMaker    *" << std::endl
+	    << "*           Start        *" << std::endl
+	    << "**************************" << std::endl << endl;
+
+  //Constants, cut and initial parameters
+  Float_t cVtxZ[2] = {-70., 70.};
+  Float_t cVtxR[2] = {0., 3.};
+  Float_t cVtxXShift = 0.;
+  Float_t cVtxYShift = 0.;
+  Float_t cVtxVpdVzDiff[2] = {-7., 7.};
+  Float_t cParticleMom[2] = {0.15, 1.65};
+  Float_t cTrackDcaGlobal[2] = {-0.01, 3.};
+  Int_t   cTrackNHits[2] = {15, 50};
+  Int_t   cTrackNHitsFit[2] = {5, 50};
+  Float_t cTrackEta[2] = {-1., 1.};
+  Int_t   cTrackFlag[2] = {0, 1000};
+  Bool_t  cZdcRefmultCorrection = true;   //Should be false for low energies
+
+  //Inclusive cuts
+  Bool_t  cSelectPions = true;
+  Bool_t  cSelectKaons = true;
+  Bool_t  cSelectTpcPions = true;
+  Bool_t  cSelectTpcKaons = true;
+  Bool_t  cSelectTofKaons = true;
+  Bool_t  cSelectTofPions = true;
+  
+
+  Float_t cPionTpcNSigmaElectron[2] = {-2., 2.};
+  Float_t cPionTpcNSigmaPion[2] = {-3., 3.};
+  Float_t cPionTpcNSigmaKaon[2] = {-2., 2.};
+  Float_t cPionTpcNSigmaProton[2] = {-3., 3.};
+  
+  Float_t cKaonTpcNSigmaElectron[2] = {-3., 3.};
+  Float_t cKaonTpcNSigmaPion[2] = {-2., 2.};
+  Float_t cKaonTpcNSigmaKaon[2] = {-3., 3.};
+  Float_t cKaonTpcNSigmaProton[2] = {-3., 3.};
+  
+  Float_t cKaonTofMom[2] = {0.55, 1.65};
+  Float_t cKaonMassSqr[2] = {0.16, 0.35};
+  Float_t cKaonTpcMom[2] = {0.15, 0.75};
+  
+  Float_t cPionTofMom[2] = {0.55, 1.65};
+  Float_t cPionMassSqr[2] = {-0.05, 0.09};
+  Float_t cPionTpcMom[2] = {0.15, 0.75};
+  
+  //Load libraries
+  std::cout << "Loading libraries..." << std::endl;
+  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
+  loadSharedLibraries();
+  gSystem->Load("libMinuit");
+  gSystem->Load("StChain");
+  gSystem->Load("StRefMultCorr");
+  gSystem->Load("StFlowMaker");
+  gSystem->Load("StHbtMaker");
+  gSystem->Load("StFemtoDstMaker");
+  gSystem->Load("libStFemtoDstMaker.so");
+  std::cout << "Libraries have been successfully loaded" << std::endl;
+
+  //Create chain
+  StChain *mChain = new StChain("StChain");
+  mChain->SetDebug(0);
+  StMuDebug::setLevel(0);
+
+  //Create MuDst reader
+  int nMaxFileToRead = 1e9;
+  std::cout << Form("Try to read file: %s", inFileList) << std::endl;
+  StMuDstMaker *muDstMaker = 
+    new StMuDstMaker(0, 0, "", inFileList, "MuDst", nMaxFileToRead);
+  muDstMaker->SetStatus("*",0);
+  muDstMaker->SetStatus("Event*",1);
+  muDstMaker->SetStatus("MuEvent*",1);
+  muDstMaker->SetStatus("PrimaryVertices*",1);
+  muDstMaker->SetStatus("PrimaryTracks*",1);
+  muDstMaker->SetStatus("GlobalTracks*",1);
+  muDstMaker->SetStatus("BTofHeader*",1);
+  TChain *mCurrentTree = muDstMaker->chain();
+  Long64_t mNEvents = mCurrentTree->GetEntries();
+  std::cout << Form("Number of events in TTree: %d", mNEvents) << std::endl;
+
+  //Heavy ion related makers
+  StRefMultCorr *mRefMult = new StRefMultCorr("refmult");
+
+  //StFemtoDstMaker initialization
+  mFemtoDstInclSelector = new StFemtoDstInclusiveSelector(muDstMaker,
+							  outFileName);
+    
+  mFemtoDstInclSelector->AddTriggerId(320000);      //Trigger selection for pp500 y2011
+  mFemtoDstInclSelector->AddTriggerId(320001);
+  mFemtoDstInclSelector->AddTriggerId(320011);
+  mFemtoDstInclSelector->AddTriggerId(320021);
+  mFemtoDstInclSelector->AddTriggerId(330021);
+  mFemtoDstInclSelector->SetVtxZCut(cVtxZ[0], cVtxZ[1]);
+  mFemtoDstInclSelector->SetVtxRCut(cVtxR[0], cVtxR[1]);
+  mFemtoDstInclSelector->SetVtxShift(cVtxXShift, cVtxYShift);
+  mFemtoDstInclSelector->SetVtxVpdVzDiffCut(cVtxVpdVzDiff[0], cVtxVpdVzDiff[1]);
+  mFemtoDstInclSelector->SetParticleMomentum(cParticleMom[0], cParticleMom[1]);
+  mFemtoDstInclSelector->SetTrackDcaGlobal(cTrackDcaGlobal[0], cTrackDcaGlobal[1]);
+  mFemtoDstInclSelector->SetTrackNHits(cTrackNHits[0], cTrackNHits[1]);
+  mFemtoDstInclSelector->SetTrackNHitsFit(cTrackNHitsFit[0], cTrackNHitsFit[1]);
+  mFemtoDstInclSelector->SetTrackEta(cTrackEta[0], cTrackEta[1]);
+  mFemtoDstInclSelector->SetTrackFlag(cTrackFlag[0], cTrackFlag[1]);
+  mFemtoDstInclSelector->SetCollisionTypeAuAu(isAuAu);
+  mFemtoDstInclSelector->SetAuAuZdcCoincidenceEnergy(cZdcRefmultCorrection);
+  mFemtoDstInclSelector->SetMuDstMaker(muDstMaker);
+
+  //Inclusive selection
+  //Pions
+  mFemtoDstInclSelector->SetSelectProtons(false);    
+  mFemtoDstInclSelector->SetSelectKaons(cSelectKaons);
+  mFemtoDstInclSelector->SetSelectPions(cSelectPions);
+  mFemtoDstInclSelector->SetSelectTpcPions(cSelectTpcPions);
+  mFemtoDstInclSelector->SetSelectTpcKaons(cSelectTpcKaons);
+  mFemtoDstInclSelector->SetSelectTofPions(cSelectTofPions);
+  mFemtoDstInclSelector->SetSelectTofKaons(cSelectTofKaons);
+  
+  mFemtoDstInclSelector->SetPionTpcMom(cPionTpcMom[0], cPionTpcMom[1]);
+  mFemtoDstInclSelector->SetPionTpcNSigmaElectron(cPionTpcNSigmaElectron[0],cPionTpcNSigmaElectron[1]);
+  mFemtoDstInclSelector->SetPionTpcNSigmaPion(cPionTpcNSigmaPion[0],cPionTpcNSigmaPion[1]);
+  mFemtoDstInclSelector->SetPionTpcNSigmaKaon(cPionTpcNSigmaKaon[0],cPionTpcNSigmaKaon[1]);
+  mFemtoDstInclSelector->SetPionTpcNSigmaProton(cPionTpcNSigmaProton[0],cPionTpcNSigmaProton[1]);
+  mFemtoDstInclSelector->SetPionTofMom(cPionTofMom[0],cPionTofMom[1]);
+  mFemtoDstInclSelector->SetPionMassSqr(cPionMassSqr[0],cPionMassSqr[1]);
+
+  mFemtoDstInclSelector->SetKaonTpcMom(cKaonTpcMom[0], cKaonTpcMom[1]);
+  mFemtoDstInclSelector->SetKaonTpcNSigmaElectron(cKaonTpcNSigmaElectron[0],cKaonTpcNSigmaElectron[1]);
+  mFemtoDstInclSelector->SetKaonTpcNSigmaPion(cKaonTpcNSigmaPion[0],cKaonTpcNSigmaPion[1]);
+  mFemtoDstInclSelector->SetKaonTpcNSigmaKaon(cKaonTpcNSigmaKaon[0],cKaonTpcNSigmaKaon[1]);
+  mFemtoDstInclSelector->SetKaonTpcNSigmaProton(cKaonTpcNSigmaProton[0],cKaonTpcNSigmaProton[1]);
+  mFemtoDstInclSelector->SetKaonTofMom(cKaonTofMom[0],cKaonTofMom[1]);
+  mFemtoDstInclSelector->SetKaonMassSqr(cKaonMassSqr[0],cKaonMassSqr[1]);
+
+  //Chain initialization and loop over events
+  mChain->Init();
+  Int_t iReturn = 0;
+  Int_t mNEventsProcessed = 0;
+  Float_t mPercentStep = 0.05;
+  Float_t mPercentCounter = mPercentStep;
+  Float_t mProgress = 0.;
+
+  time_t mStartTime, mStopTime, mDiffTime;   // time
+  float mFrac;                               // time
+  mStartTime = time(0);                      // time
+
+  for(Int_t iEvent=0; iEvent<mNEvents; iEvent++) {
+
+    mProgress = (Float_t)iEvent/(Float_t)mNEvents;
+    mNEventsProcessed++;
+    
+    if(mProgress >= mPercentCounter) {
+      mPercentCounter += mPercentStep;
+      mStopTime = time(0);
+      mDiffTime = difftime(mStopTime, mStartTime);
+      mFrac = (float)mDiffTime*(float)(mNEvents - iEvent)/(float)iEvent;
+      std::cout << Form("Processing progress: %3.0f%%. ETA: %.1f sec", mProgress*100., mFrac)
+		<< std::endl;
+    }
+
+    mChain->Clear();
+    iReturn = mChain->Make(iEvent);
+    if(iReturn != 0) {
+      std::cout << "Error has been occured. Event processing has been stopped"
+		<< std::endl;
+      break;
+    }
+  } //for(Int_t iEvent=0; iEvent<mNEvents; iEvent++)
+
+  mChain->Finish();
+  delete mFemtoDstInclSelector;
+  delete mChain;
+
+  std::cout << "**************************" << std::endl
+	    << "*       FemtoDstMaker    *" << std::endl
+	    << "*          Finish        *" << std::endl
+	    << "**************************" << std::endl << endl;
+}

+ 1 - 1
macros/femtoDst/FemtoDstMaker_pions_pp500_y2011.C

@@ -5,7 +5,7 @@
 #include "TString.h"
 #include <iostream>
 
-const Char_t* defaultInFile = "/star/data01/pwg/gnigmat/dataFiles/kk/pp500/2011/kaon_pp500_2011.MuDst.root";
+const Char_t *defaultInFile  = "/star/data01/pwg/gnigmat/dataFiles/kk/pp500/2011/kaon_pp500_2011.MuDst.root";
 const Char_t *defaultOutFile = "oFemtoDst_pp500_2011.root";
 
 class StFemtoDstInclusiveSelector;

+ 63 - 0
macros/femtoDst/FemtoDstQA_pp500_y2011.C

@@ -0,0 +1,63 @@
+#include <TFile.h>
+#include "TROOT.h"
+#include "TSystem.h"
+#include "TChain.h"
+#include "TTree.h"
+
+//_____________________
+void FemtoDstQA_pp500_y2011(const char *flist = "test.list",
+			    const char *ofile = "femtoQA_pp500_y2011.root") {
+  std::cout << "Input file list = " << flist << '\n';
+  if (!flist) {
+    std::cout << "Input file list = NULL!\n";
+    return;
+  }
+  
+  //
+  // Load libraries
+  //
+  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
+  loadSharedLibraries();
+  gSystem->Load("libMinuit");
+  gSystem->Load("StMuDSTMaker");
+  gSystem->Load("StChain");
+  gSystem->Load("StRefMultCorr");
+  gSystem->Load("StFlowMaker");
+  gSystem->Load("StFlowAnalysisMaker");
+  gSystem->Load("StHbtMaker");
+  gSystem->Load("StFemtoDstMaker");
+
+  StChain           *chain = 0;
+  StFemtoDstQAMaker *maker = 0;
+
+  //
+  // List of member links in the chain
+  //
+  chain = new StChain("StChain");
+  maker = new StFemtoDstQAMaker("", flist, ".", 12000);
+  maker->SetOutFileName(ofile);
+  maker->SetRunIdParameters(2000, 12021028, 12098009); // pp500 - production_pp500_2011
+
+  chain->Init();
+
+  //
+  // Loop over the links in the chain
+  //
+  int iret = 0;
+  int nEvents = maker->GetNEvents();
+  std::cout << "Entries in chain = " << nEvents << std::endl;
+  
+  for(int iev = 0; iev < nEvents; iev++) {
+    chain->Clear();
+    iret = chain->Make(iev);
+    if(iret) {
+      std::cout << "Error in Make occured. Error code: " << iret << std::endl;
+      std::cout << "iev = " << iev << '\n';
+      break;
+    } 
+  }
+  chain->Finish();
+  //Cleanup
+  delete maker;
+  delete chain;
+}

+ 6 - 7
macros/simulation/starHijing.C

@@ -45,15 +45,15 @@ void Hijing()
 
   // Setup collision frame, energy and beam species
   hijing->SetFrame("CMS",200.0);
-  hijing->SetBlue("Au");
-  hijing->SetYell("Au");  
+  hijing->SetBlue("p");
+  hijing->SetYell("Au");
   hijing->SetImpact(0.0, 30.0);       // Impact parameter min/max (fm)    0.   30.
 
   // Configure HIJING simulation
   HiParnt_t &hiparnt = hijing->hiparnt();
   {
     hiparnt.ihpr2(4) = 0;     // Jet quenching (1=yes/0=no)       0
-    hiparnt.ihpr2(3) = 0;     // Hard scattering (1=yes/0=no)
+    hiparnt.ihpr2(3) = 1;     // Hard scattering (1=yes/0=no)
     hiparnt.hipr1(10) = 2.0;  //    pT jet
     hiparnt.ihpr2(8)  = 10;   // Max number of jets / nucleon
     hiparnt.ihpr2(11) = 1;    // Set baryon production
@@ -65,14 +65,13 @@ void Hijing()
   // For more configuration options, see the HIJING manual
   // http://ntc0.lbl.gov/~xnwang/hijing/doc.html
 
-  primary -> AddGenerator(hijing);
-  primary -> SetCuts( 1.0E-6 , -1., -2.5, +2.5 );
-  
+  primary->AddGenerator(hijing);
+  primary->SetCuts( 0.1 , 2., -1., 1. ); // ptMin, ptMax, ymin, ymax, phiMin, phiMax, zMin, zMax
 }
 // ----------------------------------------------------------------------------
 // ----------------------------------------------------------------------------
 // ----------------------------------------------------------------------------
-void starHijing(Int_t nevents=100, Int_t rngSeed=1234, const char *fOut)
+void starHijing(Int_t nevents=1000, Int_t rngSeed=1234, const char *fOut)
 { 
   
   gROOT->ProcessLine(".L bfc.C");

+ 157 - 0
source/HistoCollector/BuildDependencies.C

@@ -0,0 +1,157 @@
+#include <TFile.h>
+#include <TGraphErrors.h>
+#include <TMultiGraph.h>
+#include <TCanvas.h>
+#include <TH1.h>
+#include <TString.h>
+#include <TMath.h>
+#include <TLegend.h>
+
+//_________________
+void BuildDependencies() {
+
+  const Char_t *mInFileName[2];
+  TFile *mInFile[2];
+  mInFileName[0] = "oKaonHBT_auau200_pid2.root";
+  mInFileName[1] = "oPionHBT_auau200_pid2.root";
+  for(Int_t i=0; i<2; i++) {
+    mInFile[i] = TFile::Open(mInFileName[i]);
+  }
+  
+  TFile *mOutFile = new TFile("oHBTdep_auau200_pid2.root","recreate");
+
+  ///////////////////////
+  //       Au+Au       //
+  ///////////////////////
+  TGraphErrors *mMultDepGrErr[2][3];   //Type, Charge
+  TGraphErrors *mMultKtGrErr[2][3][7]; //Type, Charge, Cenrality
+  TGraphErrors *mMultMtGrErr[2][3][7]; //Type, Charge, Cenrality
+  
+  //Loop over files and graphs
+  for(Int_t iFile=0; iFile<2; iFile++) {
+    for(Int_t iCharge=0; iCharge<3; iCharge++) {
+
+      mMultDepGrErr[iFile][iCharge] = (TGraphErrors*)mInFile[iFile]
+	->Get(Form("gMultDepRadiiGrErr_%d",iCharge));
+      if(iFile == 0) {
+	mMultDepGrErr[iFile][iCharge]->SetName(Form("gKaonMultDepRadiiGrErr_%d",iCharge));
+      }
+      else {
+	mMultDepGrErr[iFile][iCharge]->SetName(Form("gPionMultDepRadiiGrErr_%d",iCharge));
+      }
+
+      for(Int_t iCent=0; iCent<7; iCent++) {
+	mMultKtGrErr[iFile][iCharge][iCent] = (TGraphErrors*)mInFile[iFile]
+	  ->Get(Form("gMultKtDepRadiiGrErr_%d_%d",iCharge,(iCent+2)));
+	mMultMtGrErr[iFile][iCharge][iCent] = (TGraphErrors*)mInFile[iFile]
+	  ->Get(Form("gMultKtDepRadiiMtGrErr_%d_%d",iCharge,(iCent+2)));
+	if(iFile == 0) {
+	  mMultKtGrErr[iFile][iCharge][iCent]->
+	    SetName(Form("gKaonMultKtDepRadiiGrErr_%d_%d",iCharge,(iCent+2)));
+	  mMultMtGrErr[iFile][iCharge][iCent]->
+	    SetName(Form("gKaonMultMtDepRadiiGrErr_%d_%d",iCharge,(iCent+2)));
+	}
+	else {
+	  mMultKtGrErr[iFile][iCharge][iCent]->
+	    SetName(Form("gPionMultKtDepRadiiGrErr_%d_%d",iCharge,(iCent+2)));
+	  mMultMtGrErr[iFile][iCharge][iCent]->
+	    SetName(Form("gPionMultMtDepRadiiGrErr_%d_%d",iCharge,(iCent+2)));
+	}
+      } //for(Int_t iCent=2; iCent<9; iCent++)
+    } //for(Int_t iCharge=0; iCharge<3; iCharge++)
+  } //for(Int_t iFile=0; iFile<2; iFile++)
+
+  TMultiGraph *mgMultDep = new TMultiGraph;
+  for(Int_t iFile=0; iFile<2; iFile++) {
+    mgMultDep->Add(mMultDepGrErr[iFile][2]);
+  }
+
+  TMultiGraph *mgMultKtDep = new TMultiGraph;
+  TMultiGraph *mgMultMtDep = new TMultiGraph;
+  TMultiGraph *mgMultMtSumDep = new TMultiGraph;
+  for(Int_t iFile=0; iFile<2; iFile++) {
+
+    //Separate charges
+    for(Int_t iCharge=0; iCharge<2; iCharge++) {
+      for(Int_t iCent=0; iCent<7; iCent++) {
+	if(iFile == 0) {
+	  mgMultKtDep->Add(mMultKtGrErr[iFile][iCharge][iCent]);
+	}
+	mgMultMtDep->Add(mMultMtGrErr[iFile][iCharge][iCent]);
+      }
+    }
+
+    //Sum of charges
+    for(Int_t iCent=0; iCent<7; iCent++) {
+      mgMultMtSumDep->Add(mMultMtGrErr[iFile][2][iCent]);
+    }
+  }
+
+  Double_t mCentrality[9] = {15,31,58,98,157,234,335,430,482};
+  Double_t mCentralityErr[9] = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1};
+  Double_t mCentralityRad[9] = {1.78, 2.13, 2.45, 2.77, 3.13, 3.45, 3.78, 3.93, 4.04};
+  Double_t mCentralityRadErr[9] = {0.03, 0.04, 0.03, 0.05, 0.04, 0.04, 0.05, 0.06, 0.05};
+  TGraphErrors *mCentrSysErrTGr = new TGraphErrors(9, mCentrality, mCentralityRad, 
+						   mCentralityErr, mCentralityRadErr);
+  mCentrSysErrTGr->SetMarkerStyle(20);
+  mCentrSysErrTGr->SetMarkerColor(1);
+  mCentrSysErrTGr->SetLineWidth(2);
+  //mgMultDep->Add(mCentrSysErrTGr);
+
+  TCanvas *cCentDep = new TCanvas("cCentDep","cCentDep",800,600);
+  mgMultDep->Draw("AP");
+  mCentrSysErrTGr->Draw("[]");
+
+  TCanvas *cCentKtDep = new TCanvas("cCentKtDep","cCentKtDep",800, 600);
+  mgMultKtDep->Draw("AP");
+
+  TCanvas *cCentMtDep = new TCanvas("cCentMtDep","cCentMtDep",800,600);
+  mgMultMtDep->Draw("AP");
+
+  TCanvas *cCentMtSumDep = new TCanvas("cCentMtSumDep","cCentMtSumDep",800.600);
+  mgMultMtSumDep->Draw("APL");
+
+  /*
+  ///////////////////////
+  //         pp        //
+  ///////////////////////
+  Double_t mKaonRadiiKtVal[4] = { 1.24, 1.14, 1.12, 1.16};
+  Double_t mKaonRadiiKtValErr[4] = {0.11, 0.10, 0.11, 0.08};
+  Double_t mKaonKtVal[4] = {0.21, 0.35, 0.45, 0.66};
+  Double_t mKaonKtValErr[4] = {0.01, 0.01, 0.01, 0.01};
+  Double_t mKaonMtVal[4];
+  Double_t KAON_MASS_SQR = 0.24371698;
+  for(Int_t i=0; i<4; i++) {
+    mKaonMtVal[i] = TMath::Sqrt(KAON_MASS_SQR + mKaonKtVal[i]);
+  } 
+  TGraphErrors *mKaonMtRadii = new TGraphErrors(4, mKaonMtVal, mKaonRadiiKtVal, 
+						mKaonKtValErr, mKaonRadiiKtValErr);
+  mKaonMtRadii->SetName(Form("mKaonMtRadii"));
+  Double_t mPionRadiiKtVal[4] = {1.32, 1.26, 1.18, 1.05};
+  Double_t mPionRadiiKtValErr[4] = {0.02, 0.02, 0.02, 0.03};
+  Double_t mPionKtVal[4] = {0.25, 0.35, 0.45, 0.6};
+  Double_t mPionKtValErr[4] = {0.01, 0.01, 0.01, 0.01};
+  Double_t mPionMtVal[4];
+  Double_t PION_MASS_SQR = 0.0194797;
+  for(Int_t i=0; i<4; i++) {
+    mPionMtVal[i] = TMath::Sqrt(PION_MASS_SQR + mPionKtVal[i]);
+  } 
+  TGraphErrors *mPionMtRadii = new TGraphErrors(4, mPionMtVal, mPionRadiiKtVal, 
+						mPionKtValErr, mPionRadiiKtValErr);
+  mPionMtRadii->SetName(Form("mPionMtRadii"));
+
+  TMultiGraph *mgSmallMtDep = new TMultiGraph;
+  mgSmallMtDep->Add(mKaonMtRadii);
+  mgSmallMtDep->Add(mPionMtRadii);
+  TCanvas *cSmallMtDep = new TCanvas("cSmallMtDep","cSmallMtDep",800,600);
+  mgSmallMtDep->Draw("AP");
+  */
+
+  cCentDep->Write();
+  cCentKtDep->Write();
+  cCentMtDep->Write();
+  cCentMtSumDep->Draw();
+
+  mOutFile->Write();
+  mOutFile->Close();
+}

+ 57 - 0
source/HistoCollector/ComparePID.C

@@ -0,0 +1,57 @@
+#include <TFile.h>
+#include <TGraphErrors.h>
+#include <TMultiGraph.h>
+#include <TCanvas.h>
+#include <TH1.h>
+#include <TString.h>
+#include <TMath.h>
+#include <TLegend.h>
+
+//_________________
+void ComparePID() {
+
+  const Char_t *mInFileName[5];
+  TFile *mInFile[5];
+  for(Int_t iFile=1; iFile<5; iFile++) {
+    mInFileName[iFile] = Form("oKaonHBT_auau200_pid%d.root",iFile);
+    mInFile[iFile] = TFile::Open(mInFileName[iFile]);
+  }
+
+  Int_t mPidColor[5] = {1, 2, 3, 4, 6};
+  Double_t mMarkerSize = 1.6;
+  Int_t mMarkerStyle[5] = {20, 21, 22, 29, 34};
+  
+  TGraphErrors *mMultDepGrErr[5][3]; //[pid][charge]
+
+  TMultiGraph *mgMultDep[3];
+  TCanvas *cMultDepCanv[3];  //[charge]
+  for(Int_t iCharge=0; iCharge<3; iCharge++) {
+    mgMultDep[iCharge] = new TMultiGraph;
+  }
+
+  TFile *mOutFile = new TFile("oPIDcomparison.root");
+  for(Int_t iFile=1; iFile<5; iFile++) {
+    for(Int_t iCharge=0; iCharge<3; iCharge++) {
+
+      mMultDepGrErr[iFile][iCharge] = (TGraphErrors*)mInFile[iFile]->
+	Get(Form("gMultDepRadiiGrErr_%d",iCharge));
+      mMultDepGrErr[iFile][iCharge]->
+	SetName(Form("mMultDepGrErr_%d_%d",iFile,iCharge));
+      mMultDepGrErr[iFile][iCharge]->SetMarkerStyle(mMarkerStyle[iFile]);
+      mMultDepGrErr[iFile][iCharge]->SetMarkerSize(mMarkerSize);
+      mMultDepGrErr[iFile][iCharge]->SetLineColor(mPidColor[iFile]);
+      mMultDepGrErr[iFile][iCharge]->SetMarkerColor(mPidColor[iFile]);
+      mMultDepGrErr[iFile][iCharge]->SetLineWidth(2);
+      mgMultDep[iCharge]->Add(mMultDepGrErr[iFile][iCharge]);
+    } //for(Int_t iCharge=0; iCharge<3; iCharge++)
+  } //for(Int_t iFile=0; iFile<5; iFile++)
+
+
+  for(Int_t iCharge=0; iCharge<3; iCharge++) {
+    cMultDepCanv[iCharge] = new TCanvas(Form("cMultDepCanv%d",iCharge),
+					Form("cMultDepCanv%d",iCharge),
+					800, 600);
+    mgMultDep[iCharge]->Draw("AP");
+  } //for(Int_t iCharge=0; iCharge<3; iCharge++)
+}
+ 

+ 575 - 0
source/HistoCollector/DataProcessing.cxx

@@ -0,0 +1,575 @@
+#include "DataProcessing.h"
+#include "FemtoUtils.h"
+#include "HbtFitter.h"
+#include <TLegend.h>
+#include <TCanvas.h>
+#include <TMath.h>
+
+const Double_t MHC     = 0.197326960277;  // GeV * fm
+const Double_t MHC_SQR = 0.038937929230;  // GeV^2 * fm^2
+const Double_t INVSQRTPI = 0.564189584;   // 1/sqrt(pi)
+const Double_t PION_MASS = 0.13957;
+const Double_t PION_MASS_SQR = 0.0194797;
+const Double_t KAON_MASS = 0.493677;
+const Double_t KAON_MASS_SQR = 0.24371698;
+
+vector<double> DoFit(TH1* histo, Int_t PartType=4,
+		     Int_t SourceFunc=0, Int_t NonFemtoFunc=1,
+		     Bool_t DrawFitLegend=false, Bool_t ExcludePoints=false);
+//_________________
+DataProcessing::DataProcessing() {
+
+  mExpFileName = "~/work/star/data/expKaonHBT_auau200_2011.root";
+  mOutFileName = "~/work/star/soft/star_software/processing/auau/oKaonHBT_pp200.root";
+
+  mQinvBins = 40;
+  mQinvLow = 0.;
+  mQinvHi = 0.4;
+
+  mExpFile = NULL;
+  mOutFile = NULL;
+
+  mDetSelection = 2;
+
+  mPartType = 4;
+
+  mCorrFctnColor = 1;
+  mCorrFctnLineWidth = 2;
+  mCorrFctnMarkerStyle = 20;
+  mCorrFctnMarkerSize = 1.5;
+
+  mYaxisLo = 0.6;
+  mYaxisHi = 1.6;
+
+  mMultBins = 9;
+  mKtBins = 24;
+  mKtRanges = 6;
+
+  mSourceFunc = 0;
+  mNonFemtoFunc = 0;
+  mDrawLegend = true;
+}
+
+//_________________
+DataProcessing::~DataProcessing() {
+  /* noting to do */
+}
+
+//_________________
+void DataProcessing::SetCorrFctnStyle(Int_t color, Int_t width, 
+				      Int_t style, Double_t size) {
+  mCorrFctnColor = color;
+  mCorrFctnLineWidth = width;
+  mCorrFctnMarkerStyle = style;
+  mCorrFctnMarkerSize = size;
+}
+
+//_________________
+void DataProcessing::SetHistoStyle(TH1 *h) {
+  h->SetLineColor(mCorrFctnColor);
+  h->SetLineWidth(mCorrFctnLineWidth);
+  h->SetMarkerStyle(mCorrFctnMarkerStyle);
+  h->SetMarkerSize(mCorrFctnMarkerSize);
+  h->SetMarkerColor(mCorrFctnColor);
+  h->GetYaxis()->SetRangeUser(mYaxisLo,mYaxisHi);
+}
+
+//_________________
+void DataProcessing::OpenFiles() {
+
+  std::cout << "Opening file with experimental data: " << mExpFileName;
+  mExpFile = TFile::Open(mExpFileName);
+  std::cout << "\t [DONE]" << std::endl;
+  std::cout << "Creating output file: " << mOutFileName;
+  mOutFile = new TFile(mOutFileName,"recreate");
+  std::cout << "\t [DONE]" << std::endl;
+}
+
+//_________________
+void DataProcessing::DoProcessing() {
+
+  std::cout << "Start processing..." << std::endl;
+  OpenFiles();
+  ReadHistos();
+  ProcessHistos();
+  FitHistos();
+  MakeTGraphs();
+  Finish();
+  std::cout << "Processing is finished" << std::endl;
+}
+
+//_________________
+void DataProcessing::ReadHistos() {
+
+  std::cout << "Retrieving histograms from files...";
+
+  for(Int_t iCharge=0; iCharge<2; iCharge++) {
+    //for(Int_t iDet=0; iDet<5; iDet++) {
+    for(Int_t iMult=0; iMult<mMultBins; iMult++) {
+      for(Int_t iKt=0; iKt<mKtBins; iKt++) {
+
+	if(mPartType == 4) {
+	  //Kaons
+	  TH1F *hExpKaonKtMixNum =
+	    (TH1F*)mExpFile->Get(Form("hKaonQinvMixKt_%d_%d_%d_Num_bin_%d",
+				      iCharge,mDetSelection,iMult,iKt));
+	  TH1F *hExpKaonKtMixDen =
+	    (TH1F*)mExpFile->Get(Form("hKaonQinvMixKt_%d_%d_%d_Den_bin_%d",
+				      iCharge,mDetSelection,iMult,iKt));
+	  hQinvMixNum.push_back(hExpKaonKtMixNum);
+	  hQinvMixDen.push_back(hExpKaonKtMixDen);
+	}
+	else if(mPartType == 3) {
+	  //Pions
+	  TH1F *hExpPionKtMixNum =
+	    (TH1F*)mExpFile->Get(Form("hPionQinvMixKt_%d_%d_%d_Num_bin_%d",
+				      iCharge,mDetSelection,iMult,iKt));
+	  TH1F *hExpPionKtMixDen =
+	    (TH1F*)mExpFile->Get(Form("hPionQinvMixKt_%d_%d_%d_Den_bin_%d",
+				      iCharge,mDetSelection,iMult,iKt));
+	  hQinvMixNum.push_back(hExpPionKtMixNum);
+	  hQinvMixDen.push_back(hExpPionKtMixDen);
+	}
+
+      } //for(Int_t iKt=0; iKt<20; iKt++)
+    } //for(Int_t iMult=0; iMult<10; iMult++)
+      //} //for(Int_t iDet=0; iDet<5; iDet++)
+  } //for(Int_t iCharge=0; iCharge<2; iCharge++)
+
+  std::cout << "\t [DONE]" << std::endl;
+
+  std::cout << Form("Number of histos in numerator:   %d", hQinvMixNum.size()) 
+	    << std::endl
+	    << Form("Number of histos in denominator: %d", hQinvMixDen.size()) 
+	    << std::endl;
+}
+
+//_________________
+void DataProcessing::ProcessHistos() {
+
+  std::cout << "Processing multiplicity histograms...";
+  ProcessMultHistos();
+  std::cout << "\t [DONE]" << std::endl;
+  std::cout << "Processing Mult vs kT histograms...";
+  ProcessMultKtHistos();
+  std::cout << "\t [DONE]" << std::endl;
+}
+
+//_________________
+TH1F* DataProcessing::BuildMultCF(Int_t charge, Int_t BinLo, Int_t BinHi) {
+
+  if(charge<0 || charge>2) {
+    std::cout << "DataProcessing::BuildMultCF [WARNING] - Wrong charge: " << charge
+	      << std::endl;
+    return 0;
+  }
+  
+  if(BinLo<0 || BinHi>=mMultBins) {
+    std::cout << "DataProcessing::BuildMultCF [WARNING] - Wrong multiplicity bins: " 
+	      << BinLo << "_" << BinHi << std::endl;
+    return 0;
+  }
+
+  TH1F *hNumSum = new TH1F("hNumSum","",mQinvBins, mQinvLow, mQinvHi);
+  TH1F *hDenSum = new TH1F("hDenSum","",mQinvBins, mQinvLow, mQinvHi);
+  TH1F *hRatio;
+  TString sName = "h";
+
+  sName += "CentCF_";
+  sName += charge;
+  sName += "_";
+  sName += BinLo;
+  sName += "_";
+  sName += BinHi;
+
+  hRatio = new TH1F(sName.Data(),"",mQinvBins,mQinvLow,mQinvHi);
+
+  if(charge<2) {
+    for(Int_t iMult=BinLo; iMult<=BinHi; iMult++) {
+      for(Int_t iKt=0; iKt<mKtBins; iKt++) {
+	hNumSum->Add(hQinvMixNum[charge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
+	hDenSum->Add(hQinvMixDen[charge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
+      } //for(Int_t iKt=0; iKt<mKtBins; iKt++)
+      std::cout << std::endl;
+    } //for(Int_t iMult=BinLo; iMult<=BinHi; iMult++)
+  } //if(charge<2)
+  else {
+    for(Int_t iCharge=0; iCharge<2; iCharge++) {
+      for(Int_t iMult=BinLo; iMult<=BinHi; iMult++) {
+	for(Int_t iKt=0; iKt<mKtBins; iKt++) {
+
+	  hNumSum->Add(hQinvMixNum[iCharge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
+	  hDenSum->Add(hQinvMixDen[iCharge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
+	} //for(Int_t iKt=0; iKt<mKtBins; iKt++)
+      } //for(Int_t iMult=BinLo; iMult<=BinHi; iMult++)
+    } //for(Int_t iCharge=0; iCharge<2; iCharge++)
+  } //if(charge==2)
+
+  Make1DCorrFctn(hRatio,hNumSum,hDenSum);
+  delete hNumSum;
+  delete hDenSum;
+
+  return hRatio;
+}
+
+//_________________
+TH1F* DataProcessing::BuildMultKtCF(Int_t charge, Int_t MultBinLo, Int_t MultBinHi,
+				    Int_t KtBinLo, Int_t KtBinHi) {
+
+  if(charge<0 || charge>2) {
+    std::cout << "DataProcessing::BuildMultKtCF [WARNING] - Wrong charge: " << charge
+	      << std::endl;
+    return 0;
+  }
+  
+  if(MultBinLo<0 || MultBinHi>=mMultBins) {
+    std::cout << "DataProcessing::BuildMultKtCF [WARNING] - Wrong multiplicity bins: " 
+	      << MultBinLo << "_" << MultBinHi << std::endl;
+    return 0;
+  }
+
+  if(KtBinLo<0 || KtBinHi>=mKtBins) {
+    std::cout << "DataProcessing::BuildMultKtCF [WARNING] - Wrong kT bins: " 
+	      << KtBinLo << "_" << KtBinHi << std::endl;
+    return 0;
+  }
+
+
+  TH1F *hNumSum = new TH1F("hNumSum","",mQinvBins, mQinvLow, mQinvHi);
+  TH1F *hDenSum = new TH1F("hDenSum","",mQinvBins, mQinvLow, mQinvHi);
+  TH1F *hRatio;
+  TString sName = "h";
+
+  sName += "MultKtCF_";
+  sName += charge;
+  sName += "_";
+  sName += MultBinLo;
+  sName += "_";
+  sName += MultBinHi;
+  sName += "_";
+  sName += KtBinLo;
+  sName += "_";
+  sName += KtBinHi;
+
+  hRatio = new TH1F(sName.Data(),"",mQinvBins,mQinvLow,mQinvHi);
+
+  if(charge<2) {
+    for(Int_t iMult=MultBinLo; iMult<=MultBinHi; iMult++) {
+      for(Int_t iKt=KtBinLo; iKt<=KtBinHi; iKt++) {
+
+	hNumSum->Add(hQinvMixNum[charge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
+	hDenSum->Add(hQinvMixDen[charge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
+      } //for(Int_t iKt=KtBinLo; iKt<=KtBinHi; iKt++)
+    } //for(Int_t iMult=MultBinLo; iMult<=MultBinHi; iMult++)
+  } //if(charge<2)
+  else {
+    for(Int_t iCharge=0; iCharge<2; iCharge++) {
+      for(Int_t iMult=MultBinLo; iMult<=MultBinHi; iMult++) {
+	for(Int_t iKt=KtBinLo; iKt<=KtBinHi; iKt++) {
+
+	  hNumSum->Add(hQinvMixNum[iCharge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
+	  hDenSum->Add(hQinvMixDen[iCharge*mMultBins*mKtBins+iMult*mKtBins+iKt]);
+	} //for(Int_t iKt=KtBinLo; iKt<=KtBinHi; iKt++)
+      } //for(Int_t iMult=MultBinLo; iMult<=MultBinHi; iMult++)
+    } //for(Int_t iCharge=0; iCharge<2; iCharge++)
+  } //if(charge==2)
+
+  Make1DCorrFctn(hRatio,hNumSum,hDenSum);
+  delete hNumSum;
+  delete hDenSum;
+
+  return hRatio;
+}
+
+//_________________
+void DataProcessing::ProcessMultHistos() {
+
+  Int_t multBinLo, multBinHi;
+  for(Int_t iCharge=0; iCharge<3; iCharge++) {
+    for(Int_t iMultBin=0; iMultBin<mMultBins; iMultBin++) {
+
+      std::cout << Form("ProcessMultHistos - iCharge: %d iMultBin: %d ", 
+			iCharge, iMultBin) << std::endl;  
+      multBinLo = iMultBin; multBinHi = iMultBin;
+      hMultMixCF.push_back( BuildMultCF(iCharge, multBinLo, multBinHi) );
+      SetHistoStyle( hMultMixCF.back() );
+
+    } //for(Int_t iMultBin=0; iMultBin<mMultBins; iMultBin++)
+  } //for(Int_t iCharge=0; iCharge<3; iCharge++)
+}
+
+//_________________
+void DataProcessing::ProcessMultKtHistos() {
+
+  Int_t multBinLo, multBinHi;
+  Int_t ktBinLo, ktBinHi;
+  for(Int_t iCharge=0; iCharge<3; iCharge++) {
+    for(Int_t iMultBin=0; iMultBin<mMultBins; iMultBin++) {
+      for(Int_t iKtBin=0; iKtBin<mKtRanges; iKtBin++) {
+
+      multBinLo = iMultBin; multBinHi = iMultBin;
+      switch(iKtBin) {
+      case 0: ktBinLo = 0; ktBinHi = 5; break;
+      case 1: ktBinLo = 6; ktBinHi = 7; break;
+      case 2: ktBinLo = 8; ktBinHi = 9; break;
+	//case 3: ktBinLo = 10; ktBinHi = 18; break;
+      case 3: ktBinLo = 10; ktBinHi = 13; break;
+      case 4: ktBinLo = 14; ktBinHi = 17; break;
+      case 5: ktBinLo = 18; ktBinHi = 23; break;
+      default: ktBinLo = 0; ktBinHi = 23; break;
+      };
+      hMultKtMixCF.push_back( BuildMultKtCF(iCharge, multBinLo, multBinHi,
+					    ktBinLo, ktBinHi) );
+      SetHistoStyle( hMultKtMixCF.back() );
+
+      } //for(Int_t iKtBin=0; iKtBin<mKtRanges; iKtBin++)
+    } //for(Int_t iMultBin=0; iMultBin<mMultBins; iMultBin++)
+  } //for(Int_t iCharge=0; iCharge<3; iCharge++)
+}
+
+//_________________
+void DataProcessing::Finish() {
+
+  for(Int_t iCharge=0; iCharge<3; iCharge++) {
+    gMultDepRadiiGrErr[iCharge]->Write();
+    for(Int_t iCent=0; iCent<9; iCent++) {
+      gMultKtDepRadiiGrErr[iCharge][iCent]->Write();
+      gMultKtDepRadiiMtGrErr[iCharge][iCent]->Write();
+    }
+  }
+  mOutFile->Write();
+  mOutFile->Close();
+  mExpFile->Close();
+}
+
+//_________________
+void DataProcessing::FitHistos() {
+
+  vector<double> mFitParams;
+
+  for(unsigned int i=0; i<hMultMixCF.size(); i++) {
+    if(!mFitParams.empty()) {
+      mFitParams.clear();
+    }
+    mFitParams = DoFit( hMultMixCF[i], mPartType, mSourceFunc, 
+			mNonFemtoFunc, mDrawLegend );
+    mMultMixRadii.push_back(mFitParams.at(0));
+    mMultMixRadiiErr.push_back(mFitParams.at(1));
+    mMultMixLambda.push_back(mFitParams.at(3));
+    mMultMixLambdaErr.push_back(mFitParams.at(4));
+  } //for(unsigned int i=0; i<hMultMixCF.size(); i++)
+
+  for(unsigned int i=0; i<hMultKtMixCF.size(); i++) {
+    if(!mFitParams.empty()) {
+      mFitParams.clear();
+    }
+    mFitParams = DoFit( hMultKtMixCF[i], mPartType, mSourceFunc, 
+			mNonFemtoFunc, mDrawLegend );
+    mMultKtMixRadii.push_back(mFitParams.at(0));
+    mMultKtMixRadiiErr.push_back(mFitParams.at(1));
+    mMultKtMixLambda.push_back(mFitParams.at(3));
+    mMultKtMixLambdaErr.push_back(mFitParams.at(4));
+  } //for(unsigned int i=0; i<hMultKtMixCF.size(); i++)
+}
+
+//_________________
+void DataProcessing::MakeTGraphs() {
+
+  Double_t mMultRadiiArr[3][9];
+  Double_t mMultRadiiErrArr[3][9];
+  Double_t mCentrality[9] = {15,31,58,98,157,234,335,430,482};
+  Double_t mCentralityErr[9] = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1};
+
+  Double_t mRadiiTmp[9];
+  Double_t mRadiiErrTmp[9];
+
+  //Double_t mKtVal[4] = {0.21, 0.35, 0.45, 0.66};
+  //Double_t mKtErrVal[4] = {0.01, 0.01, 0.01, 0.01};
+  //Double_t mKtRadiiTmp[4];
+  //Double_t mKtRadiiErrTmp[4];
+  //Double_t mMtVal[4];
+  Double_t mKtRadiiTmp[6];
+  Double_t mKtRadiiErrTmp[6];
+  Double_t mKtVal[6] = {0.18, 0.31, 0.41, 0.55, 0.74, 0.97};
+  Double_t mKtErrVal[6] = {0.01, 0.01, 0.01, 0.01, 0.01, 0.01};
+  Double_t mMtVal[6];
+  if(mPartType == 3) {
+    for(Int_t i=0; i<mKtRanges; i++) {
+      mMtVal[i] = TMath::Sqrt(PION_MASS_SQR + mKtVal[i]);
+    }
+  }
+  else if(mPartType == 4) {
+    for(Int_t i=0; i<mKtRanges; i++) {
+      mMtVal[i] = TMath::Sqrt(KAON_MASS_SQR + mKtVal[i]);
+    }   
+  }
+
+  Int_t mMultDepColor[3] = {2, 4, 1};
+  Int_t mMultDepMarkerStyle;
+  Double_t mMultDepMarkerSize = 1.6;
+  if(mPartType == 4) {
+    mMultDepMarkerStyle = 20;
+  }
+  else {
+    mMultDepMarkerStyle = 21;
+  }
+
+  Int_t mMultKtDepColor[7] = {2, 4, 6, 8, 9, 28, 30};
+  Double_t mMultKtDepMarkerSize = 1.5;
+  Int_t mMultKtDepMarkerStyle[3];
+  if(mPartType == 4) {
+    mMultKtDepMarkerStyle[0] = 20;
+    mMultKtDepMarkerStyle[1] = 24;
+    mMultKtDepMarkerStyle[2] = 34;
+  }
+  else if(mPartType == 3) {
+    mMultKtDepMarkerStyle[0] = 21;
+    mMultKtDepMarkerStyle[1] = 25;
+    mMultKtDepMarkerStyle[2] = 28;    
+  }
+   
+  //Multiplicity dependence
+  for(Int_t iCharge=0; iCharge<3; iCharge++) {
+    if(iCharge==2) {
+      std::cout << std::endl;
+      std::cout << "Radii values (centrality) for sum of charges:" << std::endl;
+    }
+    for(Int_t iCent=0; iCent<mMultBins; iCent++) {
+      if(iCent==9) continue; //Skip integrated values
+      
+      mMultRadiiArr[iCharge][iCent] = mMultMixRadii.at(iCharge*mMultBins+iCent);
+      mRadiiTmp[iCent] = mMultMixRadii.at(iCharge*mMultBins+iCent);
+      mMultRadiiErrArr[iCharge][iCent] = mMultMixRadiiErr.at(iCharge*mMultBins+iCent);
+      mRadiiErrTmp[iCent] = mMultMixRadiiErr.at(iCharge*mMultBins+iCent);
+     
+      if(iCharge==2) {
+	std::cout << " " << mMultRadiiArr[iCharge][iCent];
+      }
+    } //for(Int_t iCent=0; iCent<mMultBins; iCent++)
+
+    gMultDepRadiiGrErr[iCharge] = new TGraphErrors(mMultBins,mCentrality,mRadiiTmp,
+						   mCentralityErr,mRadiiErrTmp);
+    gMultDepRadiiGrErr[iCharge]->SetName(Form("gMultDepRadiiGrErr_%d",iCharge));
+    gMultDepRadiiGrErr[iCharge]->SetMarkerColor(mMultDepColor[iCharge]);
+    gMultDepRadiiGrErr[iCharge]->SetLineColor(mMultDepColor[iCharge]);
+    gMultDepRadiiGrErr[iCharge]->SetMarkerStyle(mMultDepMarkerStyle);
+    gMultDepRadiiGrErr[iCharge]->SetMarkerSize(mMultDepMarkerSize);
+  } //for(Int_t iCharge=0; iCharge<3; iCharge++)
+
+  //Multiplicity and kT depence
+  for(Int_t iCharge=0; iCharge<3; iCharge++) {
+    if(iCharge==2) {
+      std::cout << std::endl;
+      std::cout << "Radii values (kT) for sum of charges:" << std::endl;
+    }
+    for(Int_t iCent=0; iCent<mMultBins; iCent++) {
+      if(iCharge==2) {
+	std::cout << std::endl;
+	std::cout << "Centrality bin: " << iCent << " "; 
+      }
+      if(iCent==9) continue; //Skip integrated values
+      for(Int_t iKt=0; iKt<mKtRanges; iKt++) {
+
+	mKtRadiiTmp[iKt] = mMultKtMixRadii.at(iCharge*mMultBins*mKtRanges +
+					      iCent*mKtRanges + iKt);
+	mKtRadiiErrTmp[iKt] = mMultKtMixRadiiErr.at(iCharge*mMultBins*mKtRanges +
+						    iCent*mKtRanges + iKt);
+
+	if(iCharge==2) {
+	  std::cout << " " << mKtRadiiTmp[iKt];
+	}
+      } //for(Int_t iKt=0; iKt<mKtRanges; iKt++)
+
+      gMultKtDepRadiiGrErr[iCharge][iCent] = 
+	new TGraphErrors(mKtRanges, mKtVal, mKtRadiiTmp, mKtErrVal, mKtRadiiErrTmp);
+      gMultKtDepRadiiGrErr[iCharge][iCent]->SetName(Form("gMultKtDepRadiiGrErr_%d_%d",
+							 iCharge,iCent));
+      gMultKtDepRadiiMtGrErr[iCharge][iCent] =
+	new TGraphErrors(mKtRanges, mMtVal, mKtRadiiTmp, mKtErrVal, mKtRadiiErrTmp);
+      gMultKtDepRadiiMtGrErr[iCharge][iCent]->SetName(Form("gMultKtDepRadiiMtGrErr_%d_%d",
+							   iCharge,iCent));
+      if(iCent<2) continue; //Skip first 2 point due to the mistake.
+      gMultKtDepRadiiGrErr[iCharge][iCent]->SetMarkerColor(mMultKtDepColor[iCent]);
+      gMultKtDepRadiiGrErr[iCharge][iCent]->SetLineColor(mMultKtDepColor[iCent]);
+      gMultKtDepRadiiGrErr[iCharge][iCent]->SetMarkerSize(mMultKtDepMarkerSize);
+      gMultKtDepRadiiGrErr[iCharge][iCent]->SetMarkerStyle(mMultKtDepMarkerStyle[iCharge]);
+
+      gMultKtDepRadiiMtGrErr[iCharge][iCent]->SetMarkerColor(mMultKtDepColor[iCent]);
+      gMultKtDepRadiiMtGrErr[iCharge][iCent]->SetLineColor(mMultKtDepColor[iCent]);
+      gMultKtDepRadiiMtGrErr[iCharge][iCent]->SetMarkerSize(mMultKtDepMarkerSize);
+      gMultKtDepRadiiMtGrErr[iCharge][iCent]->SetMarkerStyle(mMultKtDepMarkerStyle[iCharge]);
+    } //for(Int_t iCent=0; iCent<mMultBins; iCent++)
+  } //for(Int_t iCharge=0; iCharge<3; iCharge++)
+}
+
+//_________________
+vector<double> DoFit(TH1* histo, Int_t PartType, Int_t SourceFunc, 
+		     Int_t NonFemtoFunc, Bool_t DrawFitLegend, Bool_t mExcludePoints) {
+
+  vector<double> mRetVect;
+
+  Int_t mNbins = histo->GetXaxis()->GetNbins();
+  Double_t mFitMinVal = 0.;
+  Double_t mDeltaMaxFitVal = -0.009; // -0.005
+  Double_t mFitMaxVal = histo->GetXaxis()->GetBinUpEdge(mNbins) + mDeltaMaxFitVal;
+  Int_t mFemtoParNum = 2;       //lambda, R
+  Int_t mNonFemtoParNum = 0;
+
+  if(NonFemtoFunc == 0)
+    mNonFemtoParNum = 1;
+  else if(NonFemtoFunc == 1)
+    mNonFemtoParNum = 2;
+  else if(NonFemtoFunc >= 2 && NonFemtoFunc<=3)
+    mNonFemtoParNum = 3;
+  else if(NonFemtoFunc >= 4)
+    mNonFemtoParNum = 2;
+
+  Int_t mFullParNum = mFemtoParNum + mNonFemtoParNum;
+  HbtFitter* mHbtFitter = new HbtFitter(PartType,SourceFunc,NonFemtoFunc);
+  mHbtFitter->Set1DMaxVal(histo->GetXaxis()->GetBinUpEdge(mNbins));
+  mHbtFitter->Set1DBinsNum(mNbins);
+  TF1* mCorrFit = new TF1("mCorrFit",mHbtFitter,&HbtFitter::FitFunction,
+                          mFitMinVal,mFitMaxVal,mFullParNum,"HbtFitter","FitFunction");
+  mCorrFit->SetParameter(0, 0.5);
+  mCorrFit->SetParLimits(0, 0., 1.05);
+  mCorrFit->SetParameter(1, 15.);
+  mCorrFit->SetParLimits(1, 0., 50.);
+  mCorrFit->SetParameter(2, 1.);
+  //mCorrFit->SetParLimits(2, 0.9, 1.1);
+  mCorrFit->SetLineWidth(3);
+  mCorrFit->SetLineColor(kRed);
+  mCorrFit->SetNpx(40);
+  histo->Fit("mCorrFit","MRQE","ep");
+  mCorrFit->Draw("same");
+
+  if(SourceFunc==0) {
+    mRetVect.push_back(mCorrFit->GetParameter(1)*MHC);                 //0 R
+    mRetVect.push_back(mCorrFit->GetParError(1)*MHC);                  //1 dR
+  }
+  else if(SourceFunc==1) {
+    mRetVect.push_back(mCorrFit->GetParameter(1)*MHC*INVSQRTPI);       //0 R
+    mRetVect.push_back(mCorrFit->GetParError(1)*MHC*INVSQRTPI);        //1 dR
+  }
+  mRetVect.push_back(mCorrFit->GetParameter(0));                       //2 Lam
+  mRetVect.push_back(mCorrFit->GetParError(0));                        //3 dLam
+  mRetVect.push_back(mCorrFit->GetChisquare()/mCorrFit->GetNDF());     //4 chi2/ndf
+  mRetVect.push_back(mCorrFit->GetParameter(2));                       //5 Norm
+  mRetVect.push_back(mCorrFit->GetParError(2));                        //6 dNorm
+
+  if(DrawFitLegend) {
+    TLegend* mLegend = new TLegend(0.4,0.5,0.9,0.9);
+    mLegend->SetFillColor(kWhite);
+    mLegend->SetHeader("Fit parameters:");
+    mLegend->AddEntry("mCorrFit","Bowler-Sinyukov fit","l");
+    mLegend->AddEntry((TObject*)0, Form( "#chi^{2}/n.d.f. = %4.3f / %2d",
+					 mCorrFit->GetChisquare(),mCorrFit->GetNDF() ),"");
+    mLegend->AddEntry((TObject*)0, Form("N = %4.3f \\pm %4.3f",
+					mRetVect.at(5),mRetVect.at(6) ), "");
+    mLegend->AddEntry((TObject*)0, Form("#lambda = %4.2f \\pm %4.2f",
+					mRetVect.at(2),mRetVect.at(3) ),"");
+    mLegend->AddEntry((TObject*)0, Form("R = %4.2f \\pm %4.2f",
+					mRetVect.at(0),mRetVect.at(1)), "");
+    mLegend->Draw();
+  }
+  return mRetVect;
+}

+ 140 - 0
source/HistoCollector/DataProcessing.h

@@ -0,0 +1,140 @@
+#ifndef DataProcessing_h
+#define DataProcessing_h
+
+#include <TFile.h>
+#include <TSystem.h>
+#include <TH1.h>
+#include <TString.h>
+#include <TAttMarker.h>
+#include <TGraphErrors.h>
+#include <iostream>
+#include <vector>
+
+using namespace std;
+
+//_________________
+class DataProcessing {
+
+ public:
+
+  DataProcessing();
+  ~DataProcessing();
+
+  void SetExpFile(const Char_t *name) { mExpFileName = name; }
+  void SetOutFile(const Char_t *name) { mOutFileName = name; }
+  void DoProcessing();
+  void SetQinvBins(Int_t    qInvBins,
+		   Double_t qInvLow,
+		   Double_t qInvHi) {
+    mQinvBins = qInvBins;
+    mQinvLow = qInvLow;
+    mQinvHi = qInvHi;
+  }
+  void SetDetectorSelection(Int_t sel) { mDetSelection = sel; }
+  void SetMultiplicityBins(Int_t bins) { mMultBins = bins; }
+  void SetCorrFctnStyle(Int_t color, Int_t width, Int_t style, Double_t size);
+  void SetParticleType(Int_t type) { mPartType = type; } // 3-pion, 4-kaon
+  void SetSourceFunction(Int_t source) { mSourceFunc = source; }
+  void SetNonFemtoFunction(Int_t func) { mNonFemtoFunc = func; }
+  void SetDrawLegend(Bool_t draw) { mDrawLegend = draw; }
+
+ private:
+
+  //
+  // Functions
+  //
+  void OpenFiles();
+  void ReadHistos();
+  void ProcessHistos();
+  void ProcessMultHistos();
+  void ProcessMultKtHistos();
+  void FitHistos();
+  void MakeTGraphs();
+  TH1F *BuildMultCF(Int_t charge, Int_t BinLo, Int_t BinHi);
+  TH1F *BuildMultKtCF(Int_t charge, Int_t MultBinLo, Int_t MultBinHi,
+		      Int_t KtBinLo, Int_t KtBinHi);
+  
+  void SetHistoStyle(TH1* h);
+  void Finish();
+
+  //
+  // Variables
+  //
+  const Char_t *mExpFileName;
+  const Char_t *mOutFileName;
+
+  Int_t mQinvBins;
+  Double_t mQinvLow;
+  Double_t mQinvHi;
+
+  Int_t mPartType;
+
+  Int_t mDetSelection;   // 0-TPC, 1-TOF, 2-TPC+TOF, 3-TOF||TPC, 4-TPC+TOC(p>pthresh)||TPC(p<pthresh)
+
+  Int_t mMultBins;
+  Int_t mKtBins;
+  Int_t mKtRanges;
+
+  Int_t mCorrFctnColor;
+  Int_t mCorrFctnLineWidth;
+  Int_t mCorrFctnMarkerStyle;
+  Double_t mCorrFctnMarkerSize;
+
+  Double_t mYaxisLo;
+  Double_t mYaxisHi;
+
+  //
+  // Fits
+  //
+  Int_t mSourceFunc;
+  Int_t mNonFemtoFunc;
+  Bool_t mDrawLegend;
+
+  //
+  // Files
+  //
+  TFile *mExpFile;
+  TFile *mOutFile;
+
+  //
+  // Input histos
+  //
+  std::vector<TH1F *> hQinvMixNum;
+  std::vector<TH1F *> hQinvMixDen;
+
+  //
+  // Multiplicity dependece
+  //
+  std::vector<TH1F *> hMultMixCF;
+
+  //
+  // Mulitplicity and kT dependence
+  //
+  std::vector<TH1F *> hMultKtMixCF;
+
+  //
+  // Fitted values
+  // Multiplicity dependence
+  //
+  std::vector<double> mMultMixRadii;
+  std::vector<double> mMultMixRadiiErr;
+  std::vector<double> mMultMixLambda;
+  std::vector<double> mMultMixLambdaErr;
+
+  //
+  // Multiplicity and kT
+  //
+  std::vector<double> mMultKtMixRadii;
+  std::vector<double> mMultKtMixRadiiErr;
+  std::vector<double> mMultKtMixLambda;
+  std::vector<double> mMultKtMixLambdaErr;
+
+  //
+  // TGraphErrors
+  //
+  TGraphErrors *gMultDepRadiiGrErr[3];
+  TGraphErrors *gMultKtDepRadiiGrErr[3][9];
+  TGraphErrors *gMultKtDepRadiiMtGrErr[3][9];
+};
+
+#endif

+ 111 - 0
source/HistoCollector/FemtoUtils.h

@@ -0,0 +1,111 @@
+#ifndef FemtoUtils_h
+#define FemtoUtils_h
+
+#include <TH1.h>
+#include <TH3.h>
+
+//_________________
+void NormalizeQinv(TH1 *h) {
+
+  Int_t mNbins = h->GetNbinsX();
+  Int_t mHNbins = mNbins/2;
+
+  Double_t scale = h->Integral(1, mNbins);
+  scale -= h->Integral(1, mHNbins);
+  h->Scale(1./scale);
+}
+
+//_________________
+void Normalize1DCorrFctn(TH1 *h) {
+
+  Int_t mNbins = h->GetNbinsX();
+  Int_t mHNbins = mNbins/2;
+
+  Double_t scale = h->Integral(1, mNbins);
+  scale -= h->Integral(1, mHNbins);
+  h->Scale((Float_t)(mNbins-mHNbins)/scale);
+}
+
+//_________________
+void NormalizeQosl(TH3 *h) {
+
+  Bool_t mIsSymmetric = false;
+  Int_t mNbinsX = h->GetNbinsX();
+  Int_t mNbinsY = h->GetNbinsY();
+  Int_t mNbinsZ = h->GetNbinsZ();
+
+  if( (h->GetXaxis()->GetBinLowEdge(1) < 0.) &&
+      (h->GetYaxis()->GetBinLowEdge(1) < 0.) &&
+      (h->GetZaxis()->GetBinLowEdge(1) < 0.) ) {
+    mIsSymmetric = true;
+  }
+
+  Double_t scale = 0.;
+
+  for(Int_t iBinX=1; iBinX<mNbinsX; iBinX++) {
+    for(Int_t iBinY=1; iBinY<mNbinsY; iBinY++) {
+      for(Int_t iBinZ=1; iBinZ<mNbinsZ; iBinZ++) {
+	scale += h->GetBinContent(iBinX, iBinY, iBinZ);
+      } //for(Int_t iBinZ=1; iBinZ<mNbinsZ; iBinZ++)
+    } //for(Int_t iBinY=1; iBinY<mNbinsY; iBinY++)
+  } //for(Int_t iBinX=1; iBinX<mNbinsX; iBinX++)
+  
+  if(mIsSymmetric) {
+
+    Int_t mQNbinsX = mNbinsX/4;
+    Int_t mQNbinsY = mNbinsY/4;
+    Int_t mQNbinsZ = mNbinsZ/4;
+
+    for(Int_t iBinX=mQNbinsX; iBinX<(mNbinsX-mQNbinsX); iBinX++) {
+      for(Int_t iBinY=mQNbinsY; iBinY<(mNbinsY-mQNbinsY); iBinY++) {
+	for(Int_t iBinZ=mQNbinsZ; iBinZ<(mNbinsZ-mQNbinsZ); iBinZ++) {
+	  scale -= h->GetBinContent(iBinX, iBinY, iBinZ);
+	} //for(Int_t iBinZ=mQNbinsZ; iBinZ<(mNbinsZ-mQNbinsZ); iBinZ++)
+      } //for(Int_t iBinY=mQNbinsY; iBinY<(mNbinsY-mQNbinsY); iBinY++)
+    } //for(Int_t iBinX=mQNbinsX; iBinX<(mNbinsX-mQNbinsX); iBinX++)
+  }
+  else {
+
+    Int_t mHNbinsX = mNbinsX/2;
+    Int_t mHNbinsY = mNbinsY/2;
+    Int_t mHNbinsZ = mNbinsZ/2;
+    for(Int_t iBinX=1; iBinX<mHNbinsX; iBinX++) {
+      for(Int_t iBinY=1; iBinY<mHNbinsY; iBinY++) {
+	for(Int_t iBinZ=1; iBinZ<mHNbinsZ; iBinZ++) {
+	  scale -= h->GetBinContent(iBinX, iBinY, iBinZ);
+	} //for(Int_t iBinZ=1; iBinZ<mNbinsZ; iBinZ++)
+      } //for(Int_t iBinY=1; iBinY<mNbinsY; iBinY++)
+    } //for(Int_t iBinX=1; iBinX<mNbinsX; iBinX++)    
+  }
+
+  h->Scale(1./scale);
+}
+
+//_________________
+void Make1DCorrFctn(TH1* hCorr, TH1* hNum, TH1* hDen) {
+
+  hNum->Sumw2();
+  hDen->Sumw2();
+  NormalizeQinv(hNum);
+  NormalizeQinv(hDen);
+  hCorr->Divide(hNum,hDen);
+}
+
+//_________________
+void Correct1DCorrFctn(TH1* hCorr, TH1* hNum, TH1* hDen) {
+  
+  Normalize1DCorrFctn(hNum);
+  Normalize1DCorrFctn(hDen);
+  hCorr->Divide(hNum,hDen);
+}
+
+//_________________
+void Make3DCorrFctn(TH3* hCorr, TH3* hNum, TH3* hDen) {
+
+  hNum->Sumw2();
+  hDen->Sumw2();
+  NormalizeQosl(hNum);
+  NormalizeQosl(hDen);
+  hCorr->Divide(hNum, hDen);
+}
+#endif

File diff suppressed because it is too large
+ 1643 - 0
source/HistoCollector/HbtFitter.cxx


+ 496 - 0
source/HistoCollector/HbtFitter.h

@@ -0,0 +1,496 @@
+/*
+  This class implements the fitting basics for different
+  HBT studies. Right now it processes the 1-dimentional
+  case only. It should be easy to expand it to 3D case.
+  For the example of the implementation see TryNewFitter.C
+*/
+#ifndef HbtFitter_h
+#define HbtFitter_h
+
+#include <TH1.h>
+#include <TH2.h>
+#include <TH3.h>
+#include <TF1.h>
+#include <TF2.h>
+#include <TF3.h>
+#include <TMinuit.h>
+#include <vector>
+
+//#define MESSAGING
+//using namespace std;
+
+
+//1fm
+double KK_LONG[50] = { 0.273707, 0.697738, 0.823735, 0.881777, 0.914707,
+                       0.936038, 0.950835, 0.961148, 0.968602, 0.974788,
+                       0.979086, 0.982599, 0.985363, 0.987413, 0.989307,
+                       0.990824, 0.991961, 0.992976, 0.993674, 0.994269,
+                       0.995043, 0.995395, 0.995902, 0.996183, 0.996628,
+                       0.996828, 0.997042, 0.997301, 0.997473, 0.997612,
+                       0.997758, 0.99783,  0.998072, 0.998223, 0.998325,
+                       0.998349, 0.998441, 0.998541, 0.998638, 0.998712,
+                       0.998765, 0.998821, 0.998893, 0.998901, 0.998976,
+                       0.999042, 0.999062, 0.999104, 0.999157, 0.999139  };
+
+
+/*
+//0.5fm
+double KK_LONG[50] = { 0.268129, 0.683624, 0.807178, 0.864549, 0.8975,
+                       0.919023, 0.934211, 0.945184, 0.953503, 0.960384,
+                       0.965604, 0.97009,  0.973684, 0.976676, 0.979422,
+                       0.981683, 0.98353,  0.985301, 0.986679, 0.987884,
+                       0.989182, 0.990022, 0.990995, 0.991648, 0.992426,
+                       0.993047, 0.993507, 0.994047, 0.994508, 0.994853,
+                       0.995177, 0.995451, 0.995844, 0.996168, 0.996398,
+                       0.99654,  0.996779, 0.996912, 0.997158, 0.997283,
+                       0.997389, 0.99755,  0.997672, 0.997731, 0.997861,
+                       0.997994, 0.998054, 0.998162, 0.99824,  0.998291  };
+*/
+
+/*
+//0.75fm
+double KK_LONG[50] = { 0.270902, 0.690659, 0.815473, 0.873245, 0.906272,
+                       0.927801, 0.94292,  0.953665, 0.961657, 0.968311,
+                       0.973152, 0.977253, 0.980482, 0.983045, 0.985411,
+                       0.987326, 0.988796, 0.990185, 0.991195, 0.992056,
+                       0.993069, 0.993616, 0.994309, 0.994725, 0.995284,
+                       0.995649, 0.995945, 0.996293, 0.996556, 0.996744,
+                       0.996967, 0.997099, 0.997373, 0.997583, 0.997725,
+                       0.997793, 0.997926, 0.998039, 0.998179, 0.998258,
+                       0.998328, 0.998409, 0.998512, 0.998522, 0.998618,
+                       0.998707, 0.998736, 0.998797, 0.998855, 0.998864  };
+*/
+
+/*
+//1.25fm
+double KK_LONG[50] = { 0.276544, 0.704849, 0.831921, 0.890064, 0.922686,
+                       0.943586, 0.957802, 0.967498, 0.974251, 0.9798,
+                       0.983482, 0.98635,  0.988659, 0.990242, 0.991732,
+                       0.992899, 0.993813, 0.994582, 0.995075, 0.995513,
+                       0.996124, 0.996359, 0.996782, 0.996986, 0.997362,
+                       0.997496, 0.997649, 0.99786,  0.998001, 0.998105,
+                       0.998217, 0.99826,  0.998471, 0.998586, 0.998666,
+                       0.99867,  0.99875,  0.998818, 0.998909, 0.998984,
+                       0.999002, 0.999053, 0.999111, 0.999125, 0.999173,
+                       0.999238, 0.999252, 0.999281, 0.999324, 0.999315  };
+*/
+
+/*
+//1.5fm
+double KK_LONG[50] = { 0.279412, 0.711976, 0.839987, 0.898032, 0.930118,
+                       0.950357, 0.963763, 0.97271,  0.97868,  0.983532,
+                       0.986616, 0.988915, 0.990848, 0.992065, 0.993266,
+                       0.994181, 0.994944, 0.995577, 0.995949, 0.996313,
+                       0.996802, 0.996958, 0.997347, 0.997503, 0.997826,
+                       0.997938, 0.998043, 0.998226, 0.998341, 0.998422,
+                       0.998532, 0.998528, 0.998727, 0.998826, 0.998888,
+                       0.998898, 0.998957, 0.999003, 0.999094, 0.999157,
+                       0.999155, 0.999202, 0.999251, 0.999276, 0.999304,
+                       0.999366, 0.999384, 0.999401, 0.999435, 0.999436  };
+*/
+
+//1fm
+/*
+double PI_LONG[50] = { 0.722252, 0.907553, 0.948263, 0.965878, 0.9756,
+                       0.981806, 0.986068, 0.989019, 0.991139, 0.992895,