Nikita Ermakov лет назад: 7
Родитель
Сommit
1551e980e0
4 измененных файлов с 1009 добавлено и 0 удалено
  1. 466 0
      macros/expPionHBT_FemtoDST_1D.C
  2. 486 0
      macros/expPionHBT_MuDST_1D.C
  3. 10 0
      macros/gen_root.C
  4. 47 0
      macros/isZombie.C

+ 466 - 0
macros/expPionHBT_FemtoDST_1D.C

@@ -0,0 +1,466 @@
+#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 = 4;                //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 int    mCentBins = 9;
+const Bool_t mUseZdcCorrection = true;
+
+//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 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 = 4;
+
+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] = {-1e9., 1e9};
+const double mnSigmaKaon[2] = {-2., 2.};
+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 mTnTNSigmaElectron[2] = {-1e9, 1e9};
+const double mTnTNSigmaPion[2] = {-1e9., 1e9};
+const double mTnTNSigmaKaon[2] = {-3., 3.};
+const double mTnTNSigmaProton[2] = {-1e9, 1e9};
+
+const double mTTTMomThresh = 0.45; // 0-TPC, 1-ToF, 2-TPC+ToF, 3-ToF||TPC, 4-TPC&&TOF(P>Pthr)||TPC(P<Pthr)
+
+//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., 0.8};         //Pair relative momentum
+const double mAveSep[2] = {5., 5000.};         //Average separation of tracks within TPC
+
+//Histograms values
+const int mPairKtBins = 24;
+const double mPairKtVal[2] = {0.05, 1.25};
+const int mQinvBins = 40;                      //10 MeV/c^2 bin width
+const double mQinvVal[2] = {0., 0.4};          //Redefine for each analysis
+
+const Char_t *defaultInFile = "oFemtoDstInclusiveSelectorTest.root";
+const Char_t *defaultOutFile = "oExpPion_FemtoDST_1D.root";
+
+//_________________
+void expPionHBT_FemtoDST_1D(const Char_t *inFileList = defaultInFile,
+			    const Char_t *oFileName = defaultOutFile) {
+
+  std::cout << "**********************************" << std::endl
+	    << "*     expPionHBT_FemtoDST_1D     *" << std::endl
+	    << "*              Start             *" << std::endl
+	    << "**********************************" << std::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("libgsl");
+  gSystem->Load("libgslcblas");
+  gSystem->Load("StFemtoDstMaker");
+
+  std::cout << "Libraries have been successfully loaded" << std::endl;
+
+  //Create chain
+  StChain* mChain = new StChain("StChain");
+  mChain->SetDebug(0);
+  StMuDebug::setLevel(0);
+
+  //Set StHbtMaker and StHbtManager
+  StHbtMaker *mHbtMaker = new StHbtMaker();
+  mHbtMaker->SetDebug(0);
+  StHbtManager *mHbtManager = mHbtMaker->HbtManager();
+
+  //Set StHbtReader
+  StHbtFemtoDstReader *femtoReader = new StHbtFemtoDstReader("", inFileList, "root", -1, true);
+  femtoReader->AddTrigId(350003);          //mb-protected
+  femtoReader->AddTrigId(350013);          //mb-protected
+  femtoReader->AddTrigId(350023);          //mb-protected
+  femtoReader->AddTrigId(350033);          //mb-protected
+  femtoReader->AddTrigId(350043);          //mb-protected
+  Long64_t mNEvents = femtoReader->GetNEvents();
+  mHbtManager->SetEventReader(femtoReader);
+
+  //Set basic event cut
+  gregsEventCut *mBasicEventCut = new gregsEventCut();
+  mBasicEventCut->SetCollisionType(1);     //0-pp, 1-AuAu
+  mBasicEventCut->SetEventMult(0, 1000);
+  mBasicEventCut->SetVertZPos(mPrimVertZ[0], mPrimVertZ[1]);
+  mBasicEventCut->SetVzVpdVzDiff(mPrimVertVpdVz[0],mPrimVertVpdVz[1]);
+  mBasicEventCut->SetVertRPos(mPrimVertR);
+  mBasicEventCut->SetCentralityDefinition(mCentBins);
+  mBasicEventCut->SetEventCentralityBins(0, 8);  //For 9 centrality bins
+  mBasicEventCut->SetZdcRefMultCorrection(mUseZdcCorrection);
+
+  //Set basic kaon single-particle cut
+  gregsSimpleTrackCut *mBasicTrackCut = new gregsSimpleTrackCut();
+  mBasicTrackCut->SetCharge(mPosCharge);
+  mBasicTrackCut->SetMass(mChargedKaonMass);
+  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->SetNSigmaElectron(mnSigmaElectron[0],mnSigmaElectron[1]);
+  mBasicTrackCut->SetNSigmaPion(mnSigmaPion[0], mnSigmaPion[1]);
+  mBasicTrackCut->SetNSigmaKaon(mnSigmaKaon[0], mnSigmaKaon[1]);
+  mBasicTrackCut->SetNSigmaProton(mnSigmaProton[0], mnSigmaProton[1]);
+
+  mBasicTrackCut->SetTofMassSqr(mKaonTofMsqr[0],mKaonTofMsqr[1]);
+  mBasicTrackCut->SetTofMomentum(mTrackTofMom[0], mTrackTofMom[1]);
+
+  mBasicTrackCut->SetTnTMomentum(mTrackTnTMom[0], mTrackTnTMom[1]);
+  mBasicTrackCut->SetTnTNSigmaElectron(mTnTNSigmaElectron[0], mTnTNSigmaElectron[1]);
+  mBasicTrackCut->SetTnTNSigmaPion(mTnTNSigmaPion[0], mTnTNSigmaPion[1]);
+  mBasicTrackCut->SetTnTNSigmaKaon(mTnTNSigmaKaon[0], mTnTNSigmaKaon[1]);
+  mBasicTrackCut->SetTnTNSigmaProton(mTnTNSigmaProton[0], mTnTNSigmaProton[1]);
+
+  mBasicTrackCut->SetTTTMomentumThresh(mTTTMomThresh);
+
+  //Set basic kaon 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, mCurrentPid;
+ Int_t mCentLo, mCentHi;
+
+  ////////////////////////////////////////////
+  //                                        //
+  //       Like-sign meson analysis         //
+  //                                        //
+  ////////////////////////////////////////////
+
+  //The multiplicity analysis matrix is 2x5x10, where
+  //the first element corresponds to the particle charge: (pos, neg),
+  //the second one corresponds to the PID selection: 
+  //0-TPC, 1-ToF, 2-TPC+ToF, 3-ToF||TPC, 4-TPC&&TOF(P>Pthr)||TPC(P<Pthr),
+  //and the third one corresponds to the RefMult intervals as listed below:
+  //Centralitiy: 9 centrality bins and integrated one.
+
+  //StHbtAnalysis, cut, CF and monitor declarations
+  StHbtVertexMultAnalysis *mPionMultAnalysis[2][5][10]; 
+
+  //Cuts
+  gregsEventCut *mPionEventCut[2][5][9];
+  gregsSimpleTrackCut *mPionTrackCut[2][5][9];
+  gregsTrackPairCut *mPionPairCut[2][5][9];
+
+  //CorrFctns
+  QinvCorrFctnKt *mPionQinvMixKtCF[2][5][9];
+
+#ifdef MONITORS
+  //Cut monitors
+  fxtEventCutMonitor* mFxtPionEventCutMonitorPass[2][5][9];
+  fxtTrackCutMonitor* mFxtPionTrackCutMonitorPass[2][5][9];
+  fxtPairCutMonitor* mFxtPionPairCutMonitorPass[2][5][9];
+  fxtEventCutMonitor* mFxtPionEventCutMonitorFail[2][5][9];
+  fxtTrackCutMonitor* mFxtPionTrackCutMonitorFail[2][5][9];
+  fxtPairCutMonitor* mFxtPionPairCutMonitorFail[2][5][9];
+#endif
+
+  //Force to work with TPC+TOF pid!!! See event loop
+  for(Int_t iCharge=0; iCharge<2; iCharge++) {
+    for(Int_t iPid=2; iPid<3; iPid++) {
+      for(Int_t iMult=0; iMult<9; iMult++) {
+
+	switch(iCharge) {
+	case 0:  mLocalCharge = mPosCharge;    break;
+	case 1:  mLocalCharge = mNegCharge;    break;
+	default: mLocalCharge = mPosCharge;    break;
+	};
+
+	switch(iPid) {
+	case 0:  mCurrentPid = 0; break; //TPC only
+	case 1:  mCurrentPid = 1; break; //TOF only
+	case 2:  mCurrentPid = 2; break; //TPC+TOF
+	case 3:  mCurrentPid = 3; break; //TPC||TOF
+	case 4:  mCurrentPid = 4; break; //TPC (p<pthresh) && TOF+TPC (p>pthresh)
+	default: mCurrentPid = 4; break;
+	};
+
+        switch(iMult) {
+        case 8: mCentLo=8; mCentHi=8; mMultLow=466; mMultHi=800;  mMultBins=2; mEvents2Mix=3; break;
+        case 7: mCentLo=7; mCentHi=7; mMultLow=396; mMultHi=465;  mMultBins=2; mEvents2Mix=3; break;
+	case 6: mCentLo=6; mCentHi=6; mMultLow=281; mMultHi=395;  mMultBins=2; mEvents2Mix=3; break;
+        case 5: mCentLo=5; mCentHi=5; mMultLow=193; mMultHi=280;  mMultBins=2; mEvents2Mix=3; break;
+        case 4: mCentLo=4; mCentHi=4; mMultLow=125; mMultHi=192;  mMultBins=2;  mEvents2Mix=3; break;
+        case 3: mCentLo=3; mCentHi=3; mMultLow=76;  mMultHi=124;  mMultBins=2;  mEvents2Mix=3; break;
+	case 2: mCentLo=2; mCentHi=2; mMultLow=43;  mMultHi=75;   mMultBins=2;  mEvents2Mix=5; break;
+        case 1: mCentLo=1; mCentHi=1; mMultLow=22;  mMultHi=42;   mMultBins=2;  mEvents2Mix=10; break;
+        case 0: mCentLo=0; mCentHi=0; mMultLow=10;  mMultHi=21;   mMultBins=2;  mEvents2Mix=10; break;
+	default: mCentLo=0; mCentHi=8; mMultLow=10;  mMultHi=800;  mMultBins=10; mEvents2Mix=5; break;
+        };
+
+	/////////////////////////////////////
+	//                                 //
+	//      Like-sign pion analysis    //
+	//                                 //
+	/////////////////////////////////////
+
+	//Create analysis
+	mPionMultAnalysis[iCharge][iPid][iMult] = new StHbtVertexMultAnalysis(mPrimVertZbins,
+									      mPrimVertZ[0],
+									      mPrimVertZ[1],
+									      mMultBins,
+									      mMultLow,
+									      mMultHi);
+
+	//Create event cut
+	mPionEventCut[iCharge][iPid][iMult] = new gregsEventCut(*mBasicEventCut);
+	mPionEventCut[iCharge][iPid][iMult]->SetCentralityDefinition(mCentBins);
+	mPionEventCut[iCharge][iPid][iMult]->SetEventCentralityBins(mCentLo, mCentHi);
+	mPionEventCut[iCharge][iPid][iMult]->SetZdcRefMultCorrection(mUseZdcCorrection);
+
+#ifdef MONITORS
+	mFxtPionEventCutMonitorPass[iCharge][iPid][iMult] = 
+	  new fxtEventCutMonitor("eventPionPass", Form("_%d_%d_%d",iCharge,iPid,iMult));
+	mFxtPionEventCutMonitorFail[iCharge][iPid][iMult] = 
+	  new fxtEventCutMonitor("eventPionFail", Form("_%d_%d_%d",iCharge,iPid,iMult));
+	mPionEventCut[iCharge][iPid][iMult]->AddCutMonitor(mFxtPionEventCutMonitorPass[iCharge][iPid][iMult],
+							   mFxtPionEventCutMonitorFail[iCharge][iPid][iMult]);
+#endif
+	mPionMultAnalysis[iCharge][iPid][iMult]->SetEventCut(mPionEventCut[iCharge][iPid][iMult]);
+
+
+	//Create single-track cuts and add them to the analysis
+	mPionTrackCut[iCharge][iPid][iMult] = new gregsSimpleTrackCut(*mBasicTrackCut);
+	mPionTrackCut[iCharge][iPid][iMult]->SetMass(mChargedPionMass);
+	mPionTrackCut[iCharge][iPid][iMult]->SetCharge(mLocalCharge);
+	mPionTrackCut[iCharge][iPid][iMult]->SetDetectorSelection(mCurrentPid);
+
+	mPionTrackCut[iCharge][iPid][iMult]->SetTpcMomentum(0.15, 0.55);
+	mPionTrackCut[iCharge][iPid][iMult]->SetNSigmaElectron(-1e9,1e9);
+	mPionTrackCut[iCharge][iPid][iMult]->SetNSigmaPion(-2., 2.);
+	mPionTrackCut[iCharge][iPid][iMult]->SetNSigmaKaon(-1e9, 1e9);
+	mPionTrackCut[iCharge][iPid][iMult]->SetNSigmaProton(-1e9, 1e9);
+
+	mPionTrackCut[iCharge][iPid][iMult]->SetTofMassSqr(-0.02, 0.062);
+	mPionTrackCut[iCharge][iPid][iMult]->SetTofMomentum(0.15, 1.55);
+
+	mPionTrackCut[iCharge][iPid][iMult]->SetTnTMomentum(0.15, 1.55);
+	mPionTrackCut[iCharge][iPid][iMult]->SetTnTNSigmaElectron(-1e9, 1e9);
+	mPionTrackCut[iCharge][iPid][iMult]->SetTnTNSigmaPion(-3., 3.);
+	mPionTrackCut[iCharge][iPid][iMult]->SetTnTNSigmaKaon(-1e9, 1e9);
+	mPionTrackCut[iCharge][iPid][iMult]->SetTnTNSigmaProton(-1e9, 1e9);
+
+	mPionTrackCut[iCharge][iPid][iMult]->SetTTTMomentumThresh(0.45);
+
+#ifdef MONITORS
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult] = 
+	  new fxtTrackCutMonitor(Form("_trackPionPass_%d_%d_%d_",iCharge,iPid,iMult), mChargedPionMass);
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult] = 
+	  new fxtTrackCutMonitor(Form("_trackPionFail_%d_%d_%d_",iCharge,iPid,iMult), mChargedPionMass);
+	mPionTrackCut[iCharge][iPid][iMult]->AddCutMonitor(mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult],
+							   mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]);
+#endif
+	mPionMultAnalysis[iCharge][iPid][iMult]->SetFirstParticleCut(mPionTrackCut[iCharge][iPid][iMult]);
+	mPionMultAnalysis[iCharge][iPid][iMult]->SetSecondParticleCut(mPionTrackCut[iCharge][iPid][iMult]);
+
+	//Create pair cut
+	mPionPairCut[iCharge][iPid][iMult] = new gregsTrackPairCut(*mBasicPairCut);
+#ifdef MONITORS
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult] = 
+	  new fxtPairCutMonitor(Form("_pairPionPass_%d_%d_%d_",iCharge,iPid,iMult));
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->SetParticleMass(mChargedPionMass);
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult] = 
+	  new fxtPairCutMonitor(Form("_pairPionFail_%d_%d_%d_",iCharge,iPid,iMult));
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->SetParticleMass(mChargedPionMass);
+	mPionPairCut[iCharge][iPid][iMult]->AddCutMonitor(mFxtPionPairCutMonitorPass[iCharge][iPid][iMult],
+							  mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]);
+#endif
+	mPionMultAnalysis[iCharge][iPid][iMult]->SetPairCut(mPionPairCut[iCharge][iPid][iMult]);
+
+	//Correlation function initialization
+	mPionQinvMixKtCF[iCharge][iPid][iMult] = 
+	  new QinvCorrFctnKt(Form("hPionQinvMixKt_%d_%d_%d",iCharge,iPid,iMult),
+			     mQinvBins, mQinvVal[0], mQinvVal[1],
+			     mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+
+	mPionMultAnalysis[iCharge][iPid][iMult]->AddCorrFctn(mPionQinvMixKtCF[iCharge][iPid][iMult]);
+
+	//Set number of events to mix
+	mPionMultAnalysis[iCharge][iPid][iMult]->SetNumEventsToMix(mEvents2Mix);
+
+	//Add analysis to the manager
+	mHbtManager->AddAnalysis(mPionMultAnalysis[iCharge][iPid][iMult]);
+	
+      } //for(Int_t iMult=0; iMult<9; iMult++)
+    } //for(Int_t iPid=0; iPid<5; iPid++)
+  } //for(Int_t iCharge=0; iCharge<2; iCharge++)
+
+  //StChain initialization
+  mChain->Init();
+
+  //Start processing
+  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
+
+  //mNEvents = 15;
+  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_pp200_1D[ERROR] - iRet: " << iRet << std::endl;
+      break;
+    }
+  }
+
+  mChain->Finish();
+
+  //Create output file
+  TFile *oFile = new TFile(oFileName, "RECREATE");
+
+  //Write hbt histograms to the file
+  std::cout << "Writing histograms to the output file...";
+
+#ifdef MONITORS
+  for(Int_t iCharge=0; iCharge<2; iCharge++) {
+    for(Int_t iPid=0; iPid<5; iPid++) {
+      for(Int_t iMult=0; iMult<9; iMult++) {
+
+        //Pion analysis monitors
+        mFxtPionEventCutMonitorPass[iCharge][iPid][iMult]->VertexZ()->Write();
+	mFxtPionEventCutMonitorPass[iCharge][iPid][iMult]->VertexYvsVertexX()->Write();
+	mFxtPionEventCutMonitorPass[iCharge][iPid][iMult]->RefMult()->Write();
+	mFxtPionEventCutMonitorPass[iCharge][iPid][iMult]->VpdVzDiff()->Write();
+	mFxtPionEventCutMonitorFail[iCharge][iPid][iMult]->VertexZ()->Write();
+	mFxtPionEventCutMonitorFail[iCharge][iPid][iMult]->VertexYvsVertexX()->Write();
+	mFxtPionEventCutMonitorFail[iCharge][iPid][iMult]->RefMult()->Write();
+	mFxtPionEventCutMonitorFail[iCharge][iPid][iMult]->VpdVzDiff()->Write();
+
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->DCAGlobal()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->Nhits()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->P()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->Pt()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->NsigmaPion()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->NsigmaKaon()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->PvsDedx()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->PtvsDedx()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->Rapidity()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->PseudoRapidity()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->MassSqr()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->PvsInvBeta()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->DCAGlobal()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->Nhits()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->P()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->Pt()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->NsigmaPion()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->NsigmaKaon()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->PvsDedx()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->PtvsDedx()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->Rapidity()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->PseudoRapidity()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->MassSqr()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->PvsInvBeta()->Write();
+
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->Kt()->Write();
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->Minv()->Write();
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->Mt()->Write();
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->FMRvsQinv()->Write();
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->SLvsQinv()->Write();
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->AveSepVsQinv()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->Kt()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->Minv()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->Mt()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->FMRvsQinv()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->SLvsQinv()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->AveSepVsQinv()->Write();
+
+      } //for(Int_t iMult=0; iMult<9; iMult++)
+    } //for(Int_t iPid=0; iPid<5; iPid++)
+  } //for(Int_t iCharge=0; iCharge<2; iCharge++)
+#endif
+
+  //TPC+TOF identification only! See the loop
+  for(Int_t iCharge=0; iCharge<2; iCharge++) {
+    for(Int_t iPid=2; iPid<3; iPid++) {
+      for(Int_t iMult=0; iMult<9; iMult++) {
+	for(Int_t iKt=0; iKt<mPairKtBins; iKt++) {
+
+	  //Pion analysis
+	  mPionQinvMixKtCF[iCharge][iPid][iMult]->Numerator(iKt)->Write();
+	  mPionQinvMixKtCF[iCharge][iPid][iMult]->Denominator(iKt)->Write();
+	} //for(Int_t iKt=0; iKt<mPairKtBins; iKt++)
+      } //for(Int_t iMult=0; iMult<10; iMult++)
+    } //for(Int_t iPid=0; iPid<5; iPid++)
+  } //for(Int_t iCharge=0; iCharge<2; iCharge++
+
+  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 mChain;
+
+  std::cout << "*********************************" << std::endl
+	    << "*      expPionHBT_auau200_1D    *" << std::endl
+	    << "*              Finish           *" << std::endl
+	    << "*********************************" << std::endl << std::endl;
+}

+ 486 - 0
macros/expPionHBT_MuDST_1D.C

@@ -0,0 +1,486 @@
+#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 = 4;                //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 int    mCentBins = 9;
+const Bool_t mUseZdcCorrection = true;
+
+//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 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 = 4;
+
+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] = {-1e9., 1e9};
+const double mnSigmaKaon[2] = {-2., 2.};
+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 mTnTNSigmaElectron[2] = {-1e9, 1e9};
+const double mTnTNSigmaPion[2] = {-1e9., 1e9};
+const double mTnTNSigmaKaon[2] = {-3., 3.};
+const double mTnTNSigmaProton[2] = {-1e9, 1e9};
+
+const double mTTTMomThresh = 0.45; // 0-TPC, 1-ToF, 2-TPC+ToF, 3-ToF||TPC, 4-TPC&&TOF(P>Pthr)||TPC(P<Pthr)
+
+//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., 0.8};         //Pair relative momentum
+const double mAveSep[2] = {5., 5000.};         //Average separation of tracks within TPC
+
+//Histograms values
+const int mPairKtBins = 24;
+const double mPairKtVal[2] = {0.05, 1.25};
+const int mQinvBins = 40;                      //10 MeV/c^2 bin width
+const double mQinvVal[2] = {0., 0.4};          //Redefine for each analysis
+
+const Char_t *defaultInFile = "/star/data01/pwg/pusheax/devel/dataset/st_physics_12153012_raw_1030005.MuDst.root";
+const Char_t *defaultOutFile = "oExpPion_MuDST_1D.root";
+
+//_________________
+void expPionHBT_MuDST_1D(const Char_t *inFileList = defaultInFile,
+			   const Char_t *oFileName = defaultOutFile) {
+
+  std::cout << "*********************************" << std::endl
+	    << "*     expPionHBT_auau200_1D     *" << 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("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 muDstMaker
+  Int_t nMaxFileToRead = 1e9;
+  std::cout << Form("Try to read file: %s", inFileList) << std::endl;
+  StMuDstMaker *mMuDstMaker = 
+    new StMuDstMaker(0, 0, "", inFileList, "MuDst", nMaxFileToRead);
+
+  mMuDstMaker->SetStatus("*",0);
+  mMuDstMaker->SetStatus("Event*",1);
+  mMuDstMaker->SetStatus("MuEvent*",1);
+  mMuDstMaker->SetStatus("PrimaryVertices*",1);
+  mMuDstMaker->SetStatus("PrimaryTracks*",1);
+  mMuDstMaker->SetStatus("GlobalTracks*",1);
+  mMuDstMaker->SetStatus("BTofHeader*",1);
+
+  TChain *mCurrentTree = mMuDstMaker->chain();
+  Long64_t mNEvents = mCurrentTree->GetEntries();
+  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();
+
+  //Set StHbtReader
+  gregsStHbtMuDstReader* mHbtMuDstMakerReader = new gregsStHbtMuDstReader(mMuDstMaker);
+  mHbtMuDstMakerReader->setTrackType(1);
+  mHbtMuDstMakerReader->setReadTracks(true);
+  mHbtMuDstMakerReader->setReadV0s(false);
+  mHbtMuDstMakerReader->AddTriggerId(350003);          //mb-protected
+  mHbtMuDstMakerReader->AddTriggerId(350013);          //mb-protected
+  mHbtMuDstMakerReader->AddTriggerId(350023);          //mb-protected
+  mHbtMuDstMakerReader->AddTriggerId(350033);          //mb-protected
+  mHbtMuDstMakerReader->AddTriggerId(350043);          //mb-protected
+  mHbtMuDstMakerReader->SetProcessAllVertices(false);
+  mHbtManager->SetEventReader(mHbtMuDstMakerReader);
+
+  //Set basic event cut
+  gregsEventCut *mBasicEventCut = new gregsEventCut();
+  mBasicEventCut->SetCollisionType(1);     //0-pp, 1-AuAu
+  mBasicEventCut->SetEventMult(0, 1000);
+  mBasicEventCut->SetVertZPos(mPrimVertZ[0], mPrimVertZ[1]);
+  mBasicEventCut->SetVzVpdVzDiff(mPrimVertVpdVz[0],mPrimVertVpdVz[1]);
+  mBasicEventCut->SetVertRPos(mPrimVertR);
+  mBasicEventCut->SetCentralityDefinition(mCentBins);
+  mBasicEventCut->SetEventCentralityBins(0, 8);  //For 9 centrality bins
+  mBasicEventCut->SetZdcRefMultCorrection(mUseZdcCorrection);
+
+  //Set basic kaon single-particle cut
+  gregsSimpleTrackCut *mBasicTrackCut = new gregsSimpleTrackCut();
+  mBasicTrackCut->SetCharge(mPosCharge);
+  mBasicTrackCut->SetMass(mChargedKaonMass);
+  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->SetNSigmaElectron(mnSigmaElectron[0],mnSigmaElectron[1]);
+  mBasicTrackCut->SetNSigmaPion(mnSigmaPion[0], mnSigmaPion[1]);
+  mBasicTrackCut->SetNSigmaKaon(mnSigmaKaon[0], mnSigmaKaon[1]);
+  mBasicTrackCut->SetNSigmaProton(mnSigmaProton[0], mnSigmaProton[1]);
+
+  mBasicTrackCut->SetTofMassSqr(mKaonTofMsqr[0],mKaonTofMsqr[1]);
+  mBasicTrackCut->SetTofMomentum(mTrackTofMom[0], mTrackTofMom[1]);
+
+  mBasicTrackCut->SetTnTMomentum(mTrackTnTMom[0], mTrackTnTMom[1]);
+  mBasicTrackCut->SetTnTNSigmaElectron(mTnTNSigmaElectron[0], mTnTNSigmaElectron[1]);
+  mBasicTrackCut->SetTnTNSigmaPion(mTnTNSigmaPion[0], mTnTNSigmaPion[1]);
+  mBasicTrackCut->SetTnTNSigmaKaon(mTnTNSigmaKaon[0], mTnTNSigmaKaon[1]);
+  mBasicTrackCut->SetTnTNSigmaProton(mTnTNSigmaProton[0], mTnTNSigmaProton[1]);
+
+  mBasicTrackCut->SetTTTMomentumThresh(mTTTMomThresh);
+
+  //Set basic kaon 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, mCurrentPid;
+ Int_t mCentLo, mCentHi;
+
+  ////////////////////////////////////////////
+  //                                        //
+  //       Like-sign meson analysis         //
+  //                                        //
+  ////////////////////////////////////////////
+
+  //The multiplicity analysis matrix is 2x5x10, where
+  //the first element corresponds to the particle charge: (pos, neg),
+  //the second one corresponds to the PID selection: 
+  //0-TPC, 1-ToF, 2-TPC+ToF, 3-ToF||TPC, 4-TPC&&TOF(P>Pthr)||TPC(P<Pthr),
+  //and the third one corresponds to the RefMult intervals as listed below:
+  //Centralitiy: 9 centrality bins and integrated one.
+
+  //StHbtAnalysis, cut, CF and monitor declarations
+  StHbtVertexMultAnalysis *mPionMultAnalysis[2][5][10]; 
+
+  //Cuts
+  gregsEventCut *mPionEventCut[2][5][9];
+  gregsSimpleTrackCut *mPionTrackCut[2][5][9];
+  gregsTrackPairCut *mPionPairCut[2][5][9];
+
+  //CorrFctns
+  QinvCorrFctnKt *mPionQinvMixKtCF[2][5][9];
+
+#ifdef MONITORS
+  //Cut monitors
+  fxtEventCutMonitor* mFxtPionEventCutMonitorPass[2][5][9];
+  fxtTrackCutMonitor* mFxtPionTrackCutMonitorPass[2][5][9];
+  fxtPairCutMonitor* mFxtPionPairCutMonitorPass[2][5][9];
+  fxtEventCutMonitor* mFxtPionEventCutMonitorFail[2][5][9];
+  fxtTrackCutMonitor* mFxtPionTrackCutMonitorFail[2][5][9];
+  fxtPairCutMonitor* mFxtPionPairCutMonitorFail[2][5][9];
+#endif
+
+  //Force to work with TPC+TOF pid!!! See event loop
+  for(Int_t iCharge=0; iCharge<2; iCharge++) {
+    for(Int_t iPid=2; iPid<3; iPid++) {
+      for(Int_t iMult=0; iMult<9; iMult++) {
+
+	switch(iCharge) {
+	case 0:  mLocalCharge = mPosCharge;    break;
+	case 1:  mLocalCharge = mNegCharge;    break;
+	default: mLocalCharge = mPosCharge;    break;
+	};
+
+	switch(iPid) {
+	case 0:  mCurrentPid = 0; break; //TPC only
+	case 1:  mCurrentPid = 1; break; //TOF only
+	case 2:  mCurrentPid = 2; break; //TPC+TOF
+	case 3:  mCurrentPid = 3; break; //TPC||TOF
+	case 4:  mCurrentPid = 4; break; //TPC (p<pthresh) && TOF+TPC (p>pthresh)
+	default: mCurrentPid = 4; break;
+	};
+
+        switch(iMult) {
+        case 8: mCentLo=8; mCentHi=8; mMultLow=466; mMultHi=800;  mMultBins=2; mEvents2Mix=3; break;
+        case 7: mCentLo=7; mCentHi=7; mMultLow=396; mMultHi=465;  mMultBins=2; mEvents2Mix=3; break;
+	case 6: mCentLo=6; mCentHi=6; mMultLow=281; mMultHi=395;  mMultBins=2; mEvents2Mix=3; break;
+        case 5: mCentLo=5; mCentHi=5; mMultLow=193; mMultHi=280;  mMultBins=2; mEvents2Mix=3; break;
+        case 4: mCentLo=4; mCentHi=4; mMultLow=125; mMultHi=192;  mMultBins=2;  mEvents2Mix=3; break;
+        case 3: mCentLo=3; mCentHi=3; mMultLow=76;  mMultHi=124;  mMultBins=2;  mEvents2Mix=3; break;
+	case 2: mCentLo=2; mCentHi=2; mMultLow=43;  mMultHi=75;   mMultBins=2;  mEvents2Mix=5; break;
+        case 1: mCentLo=1; mCentHi=1; mMultLow=22;  mMultHi=42;   mMultBins=2;  mEvents2Mix=10; break;
+        case 0: mCentLo=0; mCentHi=0; mMultLow=10;  mMultHi=21;   mMultBins=2;  mEvents2Mix=10; break;
+	default: mCentLo=0; mCentHi=8; mMultLow=10;  mMultHi=800;  mMultBins=10; mEvents2Mix=5; break;
+        };
+
+	/////////////////////////////////////
+	//                                 //
+	//      Like-sign pion analysis    //
+	//                                 //
+	/////////////////////////////////////
+
+	//Create analysis
+	mPionMultAnalysis[iCharge][iPid][iMult] = new StHbtVertexMultAnalysis(mPrimVertZbins,
+									      mPrimVertZ[0],
+									      mPrimVertZ[1],
+									      mMultBins,
+									      mMultLow,
+									      mMultHi);
+
+	//Create event cut
+	mPionEventCut[iCharge][iPid][iMult] = new gregsEventCut(*mBasicEventCut);
+	mPionEventCut[iCharge][iPid][iMult]->SetCentralityDefinition(mCentBins);
+	mPionEventCut[iCharge][iPid][iMult]->SetEventCentralityBins(mCentLo, mCentHi);
+	mPionEventCut[iCharge][iPid][iMult]->SetZdcRefMultCorrection(mUseZdcCorrection);
+
+#ifdef MONITORS
+	mFxtPionEventCutMonitorPass[iCharge][iPid][iMult] = 
+	  new fxtEventCutMonitor("eventPionPass", Form("_%d_%d_%d",iCharge,iPid,iMult));
+	mFxtPionEventCutMonitorFail[iCharge][iPid][iMult] = 
+	  new fxtEventCutMonitor("eventPionFail", Form("_%d_%d_%d",iCharge,iPid,iMult));
+	mPionEventCut[iCharge][iPid][iMult]->AddCutMonitor(mFxtPionEventCutMonitorPass[iCharge][iPid][iMult],
+							   mFxtPionEventCutMonitorFail[iCharge][iPid][iMult]);
+#endif
+	mPionMultAnalysis[iCharge][iPid][iMult]->SetEventCut(mPionEventCut[iCharge][iPid][iMult]);
+
+
+	//Create single-track cuts and add them to the analysis
+	mPionTrackCut[iCharge][iPid][iMult] = new gregsSimpleTrackCut(*mBasicTrackCut);
+	mPionTrackCut[iCharge][iPid][iMult]->SetMass(mChargedPionMass);
+	mPionTrackCut[iCharge][iPid][iMult]->SetCharge(mLocalCharge);
+	mPionTrackCut[iCharge][iPid][iMult]->SetDetectorSelection(mCurrentPid);
+
+	mPionTrackCut[iCharge][iPid][iMult]->SetTpcMomentum(0.15, 0.55);
+	mPionTrackCut[iCharge][iPid][iMult]->SetNSigmaElectron(-1e9,1e9);
+	mPionTrackCut[iCharge][iPid][iMult]->SetNSigmaPion(-2., 2.);
+	mPionTrackCut[iCharge][iPid][iMult]->SetNSigmaKaon(-1e9, 1e9);
+	mPionTrackCut[iCharge][iPid][iMult]->SetNSigmaProton(-1e9, 1e9);
+
+	mPionTrackCut[iCharge][iPid][iMult]->SetTofMassSqr(-0.02, 0.062);
+	mPionTrackCut[iCharge][iPid][iMult]->SetTofMomentum(0.15, 1.55);
+
+	mPionTrackCut[iCharge][iPid][iMult]->SetTnTMomentum(0.15, 1.55);
+	mPionTrackCut[iCharge][iPid][iMult]->SetTnTNSigmaElectron(-1e9, 1e9);
+	mPionTrackCut[iCharge][iPid][iMult]->SetTnTNSigmaPion(-3., 3.);
+	mPionTrackCut[iCharge][iPid][iMult]->SetTnTNSigmaKaon(-1e9, 1e9);
+	mPionTrackCut[iCharge][iPid][iMult]->SetTnTNSigmaProton(-1e9, 1e9);
+
+	mPionTrackCut[iCharge][iPid][iMult]->SetTTTMomentumThresh(0.45);
+
+#ifdef MONITORS
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult] = 
+	  new fxtTrackCutMonitor(Form("_trackPionPass_%d_%d_%d_",iCharge,iPid,iMult), mChargedPionMass);
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult] = 
+	  new fxtTrackCutMonitor(Form("_trackPionFail_%d_%d_%d_",iCharge,iPid,iMult), mChargedPionMass);
+	mPionTrackCut[iCharge][iPid][iMult]->AddCutMonitor(mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult],
+							   mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]);
+#endif
+	mPionMultAnalysis[iCharge][iPid][iMult]->SetFirstParticleCut(mPionTrackCut[iCharge][iPid][iMult]);
+	mPionMultAnalysis[iCharge][iPid][iMult]->SetSecondParticleCut(mPionTrackCut[iCharge][iPid][iMult]);
+
+	//Create pair cut
+	mPionPairCut[iCharge][iPid][iMult] = new gregsTrackPairCut(*mBasicPairCut);
+#ifdef MONITORS
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult] = 
+	  new fxtPairCutMonitor(Form("_pairPionPass_%d_%d_%d_",iCharge,iPid,iMult));
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->SetParticleMass(mChargedPionMass);
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult] = 
+	  new fxtPairCutMonitor(Form("_pairPionFail_%d_%d_%d_",iCharge,iPid,iMult));
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->SetParticleMass(mChargedPionMass);
+	mPionPairCut[iCharge][iPid][iMult]->AddCutMonitor(mFxtPionPairCutMonitorPass[iCharge][iPid][iMult],
+							  mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]);
+#endif
+	mPionMultAnalysis[iCharge][iPid][iMult]->SetPairCut(mPionPairCut[iCharge][iPid][iMult]);
+
+	//Correlation function initialization
+	mPionQinvMixKtCF[iCharge][iPid][iMult] = 
+	  new QinvCorrFctnKt(Form("hPionQinvMixKt_%d_%d_%d",iCharge,iPid,iMult),
+			     mQinvBins, mQinvVal[0], mQinvVal[1],
+			     mPairKtBins, mPairKtVal[0], mPairKtVal[1]);
+
+	mPionMultAnalysis[iCharge][iPid][iMult]->AddCorrFctn(mPionQinvMixKtCF[iCharge][iPid][iMult]);
+
+	//Set number of events to mix
+	mPionMultAnalysis[iCharge][iPid][iMult]->SetNumEventsToMix(mEvents2Mix);
+
+	//Add analysis to the manager
+	mHbtManager->AddAnalysis(mPionMultAnalysis[iCharge][iPid][iMult]);
+	
+      } //for(Int_t iMult=0; iMult<9; iMult++)
+    } //for(Int_t iPid=0; iPid<5; iPid++)
+  } //for(Int_t iCharge=0; iCharge<2; iCharge++)
+
+  //StChain initialization
+  mChain->Init();
+
+  //Start processing
+  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
+
+  //mNEvents = 15;
+  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_pp200_1D[ERROR] - iRet: " << iRet << std::endl;
+      break;
+    }
+  }
+
+  mChain->Finish();
+
+  //Create output file
+  TFile *oFile = new TFile(oFileName, "RECREATE");
+
+  //Write hbt histograms to the file
+  std::cout << "Writing histograms to the output file...";
+
+#ifdef MONITORS
+  for(Int_t iCharge=0; iCharge<2; iCharge++) {
+    for(Int_t iPid=0; iPid<5; iPid++) {
+      for(Int_t iMult=0; iMult<9; iMult++) {
+
+        //Pion analysis monitors
+        mFxtPionEventCutMonitorPass[iCharge][iPid][iMult]->VertexZ()->Write();
+	mFxtPionEventCutMonitorPass[iCharge][iPid][iMult]->VertexYvsVertexX()->Write();
+	mFxtPionEventCutMonitorPass[iCharge][iPid][iMult]->RefMult()->Write();
+	mFxtPionEventCutMonitorPass[iCharge][iPid][iMult]->VpdVzDiff()->Write();
+	mFxtPionEventCutMonitorFail[iCharge][iPid][iMult]->VertexZ()->Write();
+	mFxtPionEventCutMonitorFail[iCharge][iPid][iMult]->VertexYvsVertexX()->Write();
+	mFxtPionEventCutMonitorFail[iCharge][iPid][iMult]->RefMult()->Write();
+	mFxtPionEventCutMonitorFail[iCharge][iPid][iMult]->VpdVzDiff()->Write();
+
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->DCAGlobal()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->Nhits()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->P()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->Pt()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->NsigmaPion()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->NsigmaKaon()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->PvsDedx()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->PtvsDedx()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->Rapidity()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->PseudoRapidity()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->MassSqr()->Write();
+	mFxtPionTrackCutMonitorPass[iCharge][iPid][iMult]->PvsInvBeta()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->DCAGlobal()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->Nhits()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->P()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->Pt()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->NsigmaPion()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->NsigmaKaon()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->PvsDedx()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->PtvsDedx()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->Rapidity()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->PseudoRapidity()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->MassSqr()->Write();
+	mFxtPionTrackCutMonitorFail[iCharge][iPid][iMult]->PvsInvBeta()->Write();
+
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->Kt()->Write();
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->Minv()->Write();
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->Mt()->Write();
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->FMRvsQinv()->Write();
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->SLvsQinv()->Write();
+	mFxtPionPairCutMonitorPass[iCharge][iPid][iMult]->AveSepVsQinv()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->Kt()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->Minv()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->Mt()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->FMRvsQinv()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->SLvsQinv()->Write();
+	mFxtPionPairCutMonitorFail[iCharge][iPid][iMult]->AveSepVsQinv()->Write();
+
+      } //for(Int_t iMult=0; iMult<9; iMult++)
+    } //for(Int_t iPid=0; iPid<5; iPid++)
+  } //for(Int_t iCharge=0; iCharge<2; iCharge++)
+#endif
+
+  //TPC+TOF identification only! See the loop
+  for(Int_t iCharge=0; iCharge<2; iCharge++) {
+    for(Int_t iPid=2; iPid<3; iPid++) {
+      for(Int_t iMult=0; iMult<9; iMult++) {
+	for(Int_t iKt=0; iKt<mPairKtBins; iKt++) {
+
+	  //Pion analysis
+	  mPionQinvMixKtCF[iCharge][iPid][iMult]->Numerator(iKt)->Write();
+	  mPionQinvMixKtCF[iCharge][iPid][iMult]->Denominator(iKt)->Write();
+	} //for(Int_t iKt=0; iKt<mPairKtBins; iKt++)
+      } //for(Int_t iMult=0; iMult<10; iMult++)
+    } //for(Int_t iPid=0; iPid<5; iPid++)
+  } //for(Int_t iCharge=0; iCharge<2; iCharge++
+
+  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 mChain;
+
+  std::cout << "*********************************" << std::endl
+	    << "*      expPionHBT_auau200_1D    *" << std::endl
+	    << "*              Finish           *" << std::endl
+	    << "*********************************" << std::endl << std::endl;
+}

+ 10 - 0
macros/gen_root.C

@@ -0,0 +1,10 @@
+void gen_root(const char *fDirOut, unsigned int counter) {
+  TRandom rnd;
+  for (unsigned int i = 0; i < counter; i++) {
+    TFile *file = new TFile(Form("%s/gen_%i.root", fDirOut, i), "CREATE");
+    TH1F *hist = new TH1F("hist", "", 10, -1., 11.);
+    hist->Fill(rnd.Rndm() * 10.);
+    hist->Write();
+    file->Close();
+  }
+}

+ 47 - 0
macros/isZombie.C

@@ -0,0 +1,47 @@
+void isZombie(const char *findPath = "./", const char *outFile = "./zombie.list") {
+  char *fName;
+
+  std::cout << "Open directory " << findPath << "... ";
+  void *pDir = gSystem->OpenDirectory(findPath);
+  if (!pDir) {
+    std::cout << "[FAIL]" << std::endl
+	      << "\texiting..." << std::endl;
+    return;
+  }
+  std::cout << "[OK]" << std::endl;
+
+
+  std::cout << "Open file " << outFile << " for writing... ";
+  FILE *fOut = fopen(outFile, "w");
+  if (!fOut) {
+    std::cout << "[FAIL]" << std::endl
+	      << "\texiting..." << std::endl;
+    return;
+  }
+  std::cout << "[OK]" << std::endl;
+  
+
+  TString fullPath("");
+  unsigned int nZombies = 0;
+  while (fName = (char *)gSystem->GetDirEntry(pDir)) {
+    if (!strstr(fName, ".root")) continue;
+
+    fullPath = findPath;
+    fullPath += '/';
+    fullPath += fName;
+    
+    TFile file(fullPath, "READ");
+    if (file.IsOpen()) {
+      std::cout << "working with file " << fName << std::endl;
+      if (file.IsZombie()) {
+	std::cout << "zombie spotted: " << fullPath << "\n";
+	nZombies++;
+	fputs(fullPath, fOut);
+      }
+    }
+  }
+
+  std::cout << "Total zombies: " << nZombies << std::endl;
+  
+  fclose(fOut);
+}