//#define GlobalTracks_type //#define DEBUG #define __TCFIT__ #define __JB_HELICES__ #define __JB_HELICES2__ //#define __OLD_DCA__ #if !defined(__CINT__) || defined(__MAKECINT__) #include "Riostream.h" #include #include "TSystem.h" #include "TMath.h" #include "TH1.h" #include "TH2.h" #include "TH3.h" #include "TProfile.h" #include "TStyle.h" #include "TF1.h" #include "TTree.h" #include "TChain.h" #include "TFile.h" #include "TNtuple.h" #include "TCanvas.h" #include "TMinuit.h" #include "TSpectrum.h" #include "TString.h" #include "TLine.h" #include "TText.h" #include "TROOT.h" #include "TList.h" #include "TPolyMarker.h" #include "StBichsel/Bichsel.h" #include "BetheBloch.h" #include "TDirIter.h" #include "TTreeIter.h" #include "TLorentzVector.h" #include "TVector3.h" #include "THelixTrack.h" #include "TRVector.h" #include "TRMatrix.h" #include "TRSymMatrix.h" #include "StDcaGeometry.h" #include "StThreeVectorF.hh" #include "StarClassLibrary/SystemOfUnits.h" #ifdef __TCFIT__ #include "StarRoot/TCFit.h" #include "StarRoot/THelixTrack.h" #endif #include "StBichsel/Bichsel.h" #else class TSystem; class TMath; class TH1; class TH2; class TH3; class TProfile; class TStyle; class TF1; class TTree; class TChain; class TFile; class TNtuple; class TCanvas; class TMinuit; class TSpectrum; class TString; class TLine; class TText; class TROOT; class TList; class TPolyMarker; class Bichsel; class BetheBloch; class TDirIter; class TTreeIter; #endif static const Double_t __PROB_SCALE__ = 1000.; static const Double_t __SIGMA_SCALE__ = 1000.; static const Double_t ame = 0.51099907e-3; static const Double_t amPi = 0.13956995; static const Double_t amP = 0.93827231; static const Double_t amK = 0.4936770; static const Double_t pTCut = 0.1; // 0.1; static const Double_t mKpiMin = 1.4; static const Double_t mKpiMax = 2.4; static const Double_t DcaCut = 0.0600; class Bichsel; Bichsel *m_Bichsel = 0; #ifdef __CuCu__ const Int_t NTriggerIDs = 42; #else const Int_t NTriggerIDs = 32; #endif //________________________________________________________________________________ Bool_t isTrigger(UInt_t id, const UInt_t *ids) { for (UInt_t i = 0; i < NTriggerIDs; i++) { if (ids[i] && id == ids[i]) return kTRUE; } return kFALSE; } //________________________________________________________________________________ //void myMuKpi(const Char_t *files="*.MuDst.root", //const Char_t *Out="MuKpi.root") { void myMuKpiPureD0sHelicesTCFit2(const Char_t *listName = "3kPureD0MuDst.list", const Char_t *Out = "MuKpi3kPureD0sHelicesLinearGlobal2.root" #ifdef __TCFIT__ , Double_t dLCut = 2 #endif ){ if (!m_Bichsel) { gSystem->Load("StBichsel"); gSystem->Load("StarClassLibrary"); m_Bichsel = Bichsel::Instance(); } //TDirIter Dir(files); //Char_t *file = 0; //Char_t *file1 = 0; //Int_t NFiles = 0; //TTreeIter iter("MuDst"); //while ((file = (Char_t *) Dir.NextFile())) {iter.AddFile(file); NFiles++; file1 = file;} //cout << files << "\twith " << NFiles << " files" << endl; //if (! file1 ) return; //TString output(Out); //if (output == "") { //TString dir = gSystem->BaseName(file1); //dir.ReplaceAll(".MuDst",""); //output += dir; //} TFile *fOut = new TFile(Out,"recreate"); Float_t tuple[50]; for(int counter=0;counter<50;counter++){ tuple[counter]=0.0; } TNtuple *nTuple = new TNtuple("MuKpiTree","MuKpiTree","FileId:EventId:RunId:PVtxX:PVtxY:PVtxZ:PtD0:PxD0:PyD0:PzD0:EtaD0:MassD0:DvtxX:DvtxY:DvtxZ:cosThetaStar:TrackK:PtK:PzK:EtaK:DcaXYK:DcaZK:SigmaXYK:SigmaZK:TrackPi:PtPi:PzPi:EtaPi:DcaXYPi:DcaZPi:SigmaXYPi:SigmaZPi:NGlob:NPrim:DeltaPxD:DeltaPyD:DeltaPzD:DeltaEtaD:DeltaPhiD:decayLengthD:DeltaDecayLengthD"); //tuple[0] = 12.0; //nTuple->Fill(tuple); Int_t NevProc = 0; Int_t currentEvent = 0; Int_t D0found = 0; FILE *readMe = fopen("READMEmuDst.txt","w");{ // create file with the layout for data stored in text file fprintf(readMe,"The data stored in each row corresponds to data for the parent and daughter particles for one event. \n"); fprintf(readMe,"Each variable will be separated by a tab space. The list is as follows: \n"); fprintf(readMe,"FileId \t EventId \t RunId \t PrimVertexX \t PrimVertexY \t PrimVertexZ \t PtD \t PxD \t PyD \t PzD \t EtaD \t MD \t DvtxX \t DvtxY \t DvtxZ \t CosThetaStar \t KeyKGlobal \t PtK \t PzK \t EtaK \t dcaXYK \t dcaZK \t KeyPiGlobal \t PtPi \t PzPi \t EtaPi dcaXYPi \t dcaZPi \t noTracksGl \t noTracksPrim \t decayLength \t sigmaXYK \t sigmaXYPi \t slength pxK \t pyK \t pxPi \t pyPi \t StopRD \t decayOrig.X() \t decayOrig.Y() \t decayOrig.Z() \t decay.x() \t decay.y() \t decay.z() \t decayLengthHel1 \t stopRDHel1 \t decayC.x() \t decayC.y() \t decayC.z() \t decayLengthHel2 \t stopRDHel2 \t decayTCFit.X() \t decayTCFit.Y() \t decayTCFit.Z() \t decayLengthTCFit \t stopRDTCFit \n"); } fclose(readMe); cout << "**************** README Written ***************" << endl; FILE *muDstData = fopen("muDstData3kPureD0sHelicesLinearTCFit2.txt","w");{ // create file with the layout for data stored in text file //fprintf(muDstData,"FileId \t EventId \t ParentId \t ParentEta \t ParentPt \t ParentPz \t GeantId1 \t Eta1 \t Pt1 \t Pz1 \t RecoKey1 \t GeantId2 \t Eta2 \t Pt2 \t Pz2 \t RecoKey2 \n",); //fprintf(muDstData,"This is a test! \n"); //cout << "Output for " << output << endl; ifstream keysFile; keysFile.open("miniMcKeys3kPureD0s.txt"); ifstream inFile; Int_t NFiles = 0; Int_t MAXFILES = 0; Int_t MAXEVENTS = 0; Int_t fileIdIntMuDst = 0; inFile.open(listName); cout << listName << endl; while(NFiles < 3){ TTreeIter iter("MuDst"); TString t; inFile >> t; cout << "NFiles = " << NFiles << " t = " << t << endl; TString fileIdMuDst(t); fileIdMuDst.ReplaceAll("/star/institutions/ksu/fisyak/D0B/rcf1272_",""); fileIdMuDst.ReplaceAll("/star/institutions/ksu/bouchet/JAIBY/simu/1-1000/test1D0_3kev_new.MuDst.root",""); fileIdMuDst.ReplaceAll("/star/institutions/ksu/bouchet/JAIBY/simu/1001-2000/test1D0_3kev_new.MuDst.root",""); fileIdMuDst.ReplaceAll("/star/institutions/ksu/bouchet/JAIBY/simu/2001-3000/test1D0_3kev_new.MuDst.root",""); fileIdMuDst.ReplaceAll("_400evts_37.MuDst.root",""); fileIdMuDst.ReplaceAll("LD0sgoodDev_100evts_",""); fileIdMuDst.ReplaceAll("MuDst.root",""); fileIdIntMuDst = atoi(fileIdMuDst); if(!inFile.good()) break; iter.Reset(); iter.AddFile(t); NFiles++; cout << "FileId = " << fileIdIntMuDst << endl; MAXFILES = NFiles; MAXEVENTS = NFiles*400; cout << "maxfiles = " << MAXFILES << endl; cout << "\n fileList:" << listName << endl; //Declare track parameters const Int_t*& EventId =iter("MuEvent.mEventInfo.mId"); const Int_t*& RunId = iter("MuEvent.mRunInfo.mRunId"); const Int_t*& uTime = iter("MuEvent.mEventInfo.mTime"); #ifdef __TCFIT__ const Double_t*&mMagneticFieldZ = iter("MuEvent.mRunInfo.mMagneticFieldZ"); #endif const Float_t*& MuEvent_mRunInfo_mCenterOfMassEnergy = iter("MuEvent.mRunInfo.mCenterOfMassEnergy"); #ifdef __OLD_DCA__ const UInt_t*& mNTriggerId_mId = iter(Form("MuEvent.mTriggerIdCollection.mNTriggerId.mId[%i]",NTriggerIDs)); #endif const Int_t &NoPrimVertex = iter("PrimaryVertices"); const Float_t *&PrimVertexX = iter("PrimaryVertices.mPosition.mX1"); const Float_t *&PrimVertexY = iter("PrimaryVertices.mPosition.mX2"); const Float_t *&PrimVertexZ = iter("PrimaryVertices.mPosition.mX3"); const Float_t *&PrimVerSigX = iter("PrimaryVertices.mPosError.mX1"); const Float_t *&PrimVerSigY = iter("PrimaryVertices.mPosError.mX2"); const Float_t *&PrimVerSigZ = iter("PrimaryVertices.mPosError.mX3"); const Int_t &NoTracks = iter("PrimaryTracks"); const Int_t *&PrimaryTracks_mVertexIndex = iter("PrimaryTracks.mVertexIndex"); const UChar_t *&mNHitsFitInner = iter("PrimaryTracks.mNHitsFitInner"); const UChar_t *&mNHitsPossInner = iter("PrimaryTracks.mNHitsPossInner"); const UChar_t *&mNHitsFit = iter("PrimaryTracks.mNHitsFit"); const Float_t *&pT = iter("PrimaryTracks.mPt"); const Float_t *&pX = iter("PrimaryTracks.mP.mX1"); const Float_t *&pY = iter("PrimaryTracks.mP.mX2"); const Float_t *&pZ = iter("PrimaryTracks.mP.mX3"); const Float_t *&Eta = iter("PrimaryTracks.mEta"); // const Float_t *&Phi = iter("PrimaryTracks.mPhi"); #ifdef __OLD_DCA__ const Float_t *&mSigmaDcaD = iter("PrimaryTracks.mSigmaDcaD"); const Float_t *&mSigmaDcaZ = iter("PrimaryTracks.mSigmaDcaZ"); #endif const Int_t *&mIndex2Global = iter("PrimaryTracks.mIndex2Global"); // const Float_t*& mDCA_mX1 = iter("PrimaryTracks.mDCA.mX1"); // const Float_t*& mDCA_mX2 = iter("PrimaryTracks.mDCA.mX2"); // const Float_t*& mDCA_mX3 = iter("PrimaryTracks.mDCA.mX3"); // const Float_t*& mDCAGlobal_mX1 = iter("PrimaryTracks.mDCAGlobal.mX1"); // const Float_t*& mDCAGlobal_mX2 = iter("PrimaryTracks.mDCAGlobal.mX2"); // const Float_t*& mDCAGlobal_mX3 = iter("PrimaryTracks.mDCAGlobal.mX3"); // Global const Int_t &NoTracksGl = iter("GlobalTracks"); const UChar_t *&mNHitsFitInnerGl = iter("GlobalTracks.mNHitsFitInner"); const UChar_t *&mNHitsPossInnerGl = iter("GlobalTracks.mNHitsPossInner"); const UChar_t *&mNHitsFitGl = iter("GlobalTracks.mNHitsFit"); const Float_t *&pTGl = iter("GlobalTracks.mPt"); const Float_t *&pXGl = iter("GlobalTracks.mP.mX1"); const Float_t *&pYGl = iter("GlobalTracks.mP.mX2"); const Float_t *&pZGl = iter("GlobalTracks.mP.mX3"); const Float_t *&EtaGl = iter("GlobalTracks.mEta"); // const Float_t *&PhiGl = iter("GlobalTracks.mPhi"); #ifdef __OLD_DCA__ const Float_t *&mSigmaDcaDGl = iter("GlobalTracks.mSigmaDcaD"); const Float_t *&mSigmaDcaZGl = iter("GlobalTracks.mSigmaDcaZ"); #else const Int_t*& GlobalTracks_mIndex2Cov = iter("GlobalTracks.mIndex2Cov"); const Int_t& NoCovGlobTrack = iter("CovGlobTrack"); const Float_t*& CovGlobTrack_mImp = iter("CovGlobTrack.mImp"); const Float_t*& CovGlobTrack_mZ = iter("CovGlobTrack.mZ"); const Float_t*& CovGlobTrack_mPsi = iter("CovGlobTrack.mPsi"); const Float_t*& CovGlobTrack_mPti = iter("CovGlobTrack.mPti"); const Float_t*& CovGlobTrack_mTan = iter("CovGlobTrack.mTan"); const Float_t*& CovGlobTrack_mCurv = iter("CovGlobTrack.mCurv"); const Float_t*& CovGlobTrack_mImpImp = iter("CovGlobTrack.mImpImp"); const Float_t*& CovGlobTrack_mZImp = iter("CovGlobTrack.mZImp"); const Float_t*& CovGlobTrack_mZZ = iter("CovGlobTrack.mZZ"); const Float_t*& CovGlobTrack_mPsiImp = iter("CovGlobTrack.mPsiImp"); const Float_t*& CovGlobTrack_mPsiZ = iter("CovGlobTrack.mPsiZ"); const Float_t*& CovGlobTrack_mPsiPsi = iter("CovGlobTrack.mPsiPsi"); const Float_t*& CovGlobTrack_mPtiImp = iter("CovGlobTrack.mPtiImp"); const Float_t*& CovGlobTrack_mPtiZ = iter("CovGlobTrack.mPtiZ"); const Float_t*& CovGlobTrack_mPtiPsi = iter("CovGlobTrack.mPtiPsi"); const Float_t*& CovGlobTrack_mPtiPti = iter("CovGlobTrack.mPtiPti"); const Float_t*& CovGlobTrack_mTanImp = iter("CovGlobTrack.mTanImp"); const Float_t*& CovGlobTrack_mTanZ = iter("CovGlobTrack.mTanZ"); const Float_t*& CovGlobTrack_mTanPsi = iter("CovGlobTrack.mTanPsi"); const Float_t*& CovGlobTrack_mTanPti = iter("CovGlobTrack.mTanPti"); const Float_t*& CovGlobTrack_mTanTan = iter("CovGlobTrack.mTanTan"); #endif // const Float_t*& mDCA_mX1Gl = iter("GlobalTracks.mDCA.mX1"); // const Float_t*& mDCA_mX2Gl = iter("GlobalTracks.mDCA.mX2"); // const Float_t*& mDCA_mX3Gl = iter("GlobalTracks.mDCA.mX3"); const Float_t*& mDCAGlobal_mX1Gl = iter("GlobalTracks.mDCAGlobal.mX1"); const Float_t*& mDCAGlobal_mX2Gl = iter("GlobalTracks.mDCAGlobal.mX2"); const Float_t*& mDCAGlobal_mX3Gl = iter("GlobalTracks.mDCAGlobal.mX3"); // const Float_t*& GlobalTracks_mFirstPoint_mX1 = iter("GlobalTracks.mFirstPoint.mX1"); // const Float_t*& GlobalTracks_mFirstPoint_mX2 = iter("GlobalTracks.mFirstPoint.mX2"); // const Float_t*& GlobalTracks_mFirstPoint_mX3 = iter("GlobalTracks.mFirstPoint.mX3"); // const Float_t*& GlobalTracks_mLastPoint_mX1 = iter("GlobalTracks.mLastPoint.mX1"); // const Float_t*& GlobalTracks_mLastPoint_mX2 = iter("GlobalTracks.mLastPoint.mX2"); // const Float_t*& GlobalTracks_mLastPoint_mX3 = iter("GlobalTracks.mLastPoint.mX3"); #ifdef __JB_HELICES__ //JB: for helices const Short_t*& QGl = iter("GlobalTracks.mHelix.mQ"); const Float_t*& OriginX = iter("GlobalTracks.mHelix.mOrigin.mX1"); const Float_t*& OriginY = iter("GlobalTracks.mHelix.mOrigin.mX2"); const Float_t*& OriginZ = iter("GlobalTracks.mHelix.mOrigin.mX3"); const Float_t*& hpX = iter("GlobalTracks.mHelix.mP.mX1"); const Float_t*& hpY = iter("GlobalTracks.mHelix.mP.mX2"); const Float_t*& hpZ = iter("GlobalTracks.mHelix.mP.mX3"); const Float_t*& hB = iter("GlobalTracks.mHelix.mB"); #endif #if 0 const Int_t *&NSigmaElectronGl = iter("GlobalTracks.mNSigmaElectron"); const Int_t *&NSigmaPionGl = iter("GlobalTracks.mNSigmaPion"); const Int_t *&NSigmaKaonGl = iter("GlobalTracks.mNSigmaKaon"); const Int_t *&NSigmaProtonGl = iter("GlobalTracks.mNSigmaProton"); #endif const Float_t*& dEdxTrackLength = iter("GlobalTracks.mProbPidTraits.mdEdxTrackLength"); const Float_t*& GdEdx = iter("GlobalTracks.mdEdx"); const Int_t NpT = 5; const Int_t Neta = 3; const Int_t Nsys = 1; const Int_t NZ = 4; const Int_t NL = 10; const Int_t NGJ = 3; const Int_t NT = 5; const Int_t NF = 13; const Int_t Ncut = 6; const Double_t pTmin[NpT] = {0.0, 0.5, 1.0, 2.0, 5.0}; const Double_t etamin[Neta] = {0.0, 0.5, 1.0}; static const Double_t L = 34.42; // half a SSD ladder's length static const Double_t R = 22.80; // Fitter Int_t NevProc = 0; #ifdef __TCFIT__ TCFitV0 dat; TCFit tc("Fit decay length",&dat); #ifndef DEBUG tc.SetDebug(0); #endif /* !DEBUG */ #endif Int_t FileId = 0; Int_t mEventId = 0; Int_t mRunId = 0; Int_t GeantK = 0; Float_t mVertexX = 0.0; Float_t mVertexY = 0.0; Float_t mVertexZ = 0.0; Float_t StopRDMc = 0.0; Float_t PtDMc = 0.0; Float_t PxDMc = 0.0; Float_t PyDMc = 0.0; Float_t PzDMc = 0.0; Float_t EtaDMc = 0.0; Float_t PhiDMc = 0.0; Float_t StopRK = 0; Int_t KeyK = 0; Int_t RecoKeyK = 0; Float_t PtKPr = 0.0; Float_t PxKPr = 0.0; Float_t PyKPr = 0.0; Float_t PzKPr = 0.0; Float_t EtaKPr = 0.0; Float_t PhiKPr = 0.0; Float_t DCAKPr = 0.0; Float_t DCAxyKPr = 0.0; Float_t DCAzKPr = 0.0; Float_t PtKGl = 0.0; Float_t PxKGl = 0.0; Float_t PyKGl = 0.0; Float_t PzKGl = 0.0; Float_t EtaKGl = 0.0; Float_t PhiKGl = 0.0; Float_t DCAKGl = 0.0; Float_t DCAxyKGl = 0.0; Float_t DCAzKGl = 0.0; Int_t GeantPi = 0; Float_t StopRPi = 0.0; Int_t KeyPi = 0; Int_t RecoKeyPi = 0; Float_t PtPiPr = 0.0; Float_t PxPiPr = 0.0; Float_t PyPiPr = 0.0; Float_t PzPiPr = 0.0; Float_t EtaPiPr = 0.0; Float_t PhiPiPr = 0.0; Float_t DCAPiPr = 0.0; Float_t DCAxyPiPr = 0.0; Float_t DCAzPiPr = 0.0; Float_t PtPiGl = 0.0; Float_t PxPiGl = 0.0; Float_t PyPiGl = 0.0; Float_t PzPiGl = 0.0; Float_t EtaPiGl = 0.0; Float_t PhiPiGl = 0.0; Float_t DCAPiGl = 0.0; Float_t DCAxyPiGl = 0.0; Float_t DCAzPiGl = 0.0; Int_t trackK = 0; Int_t trackPi = 0; keysFile >> FileId; keysFile >> mEventId; keysFile >> mRunId; keysFile >> mVertexX; keysFile >> mVertexY; keysFile >> mVertexZ; keysFile >> StopRDMc; keysFile >> PtDMc; keysFile >> PxDMc; keysFile >> PyDMc; keysFile >> PzDMc; keysFile >> EtaDMc; keysFile >> PhiDMc; keysFile >> GeantK; keysFile >> StopRK; keysFile >> KeyK; keysFile >> RecoKeyK; keysFile >> PtKPr; keysFile >> PxKPr; keysFile >> PyKPr; keysFile >> PzKPr; keysFile >> EtaKPr; keysFile >> PhiKPr; keysFile >> DCAKPr; keysFile >> DCAxyKPr; keysFile >> DCAzKPr; keysFile >> PtKGl; keysFile >> PxKGl; keysFile >> PyKGl; keysFile >> PzKGl; keysFile >> EtaKGl; keysFile >> PhiKGl; keysFile >> DCAKGl; keysFile >> DCAxyKGl; keysFile >> DCAzKGl; keysFile >> GeantPi; keysFile >> StopRPi; keysFile >> KeyPi; keysFile >> RecoKeyPi; keysFile >> PtPiPr; keysFile >> PxPiPr; keysFile >> PyPiPr; keysFile >> PzPiPr; keysFile >> EtaPiPr; keysFile >> PhiPiPr; keysFile >> DCAPiPr; keysFile >> DCAxyPiPr; keysFile >> DCAzPiPr; keysFile >> PtPiGl; keysFile >> PxPiGl; keysFile >> PyPiGl; keysFile >> PzPiGl; keysFile >> EtaPiGl; keysFile >> PhiPiGl; keysFile >> DCAPiGl; keysFile >> DCAxyPiGl; keysFile >> DCAzPiGl; Bool_t foundEvent = 0; // Fitter // loop over events Int_t nEventsMax = 1000; //Int_t nPairs = 0; //Int_t nPairsMax = 1; while(iter.Next()){ //while(nPairs < nPairsMax){ if(mEventId > EventId[0]){ currentEvent++; iter.Next(); NevProc++; } else if (mEventId < EventId[0]){ keysFile >> FileId; keysFile >> mEventId; //cout << "Find mEventId = " << mEventId << endl; keysFile >> mRunId; keysFile >> mVertexX; keysFile >> mVertexY; keysFile >> mVertexZ; keysFile >> StopRDMc; keysFile >> PtDMc; keysFile >> PxDMc; keysFile >> PyDMc; keysFile >> PzDMc; keysFile >> EtaDMc; keysFile >> PhiDMc; keysFile >> GeantK; keysFile >> StopRK; keysFile >> KeyK; keysFile >> RecoKeyK; //cout << "Find Track RecoKeyK-1 = " << RecoKeyK-1 << endl; keysFile >> PtKPr; keysFile >> PxKPr; keysFile >> PyKPr; keysFile >> PzKPr; keysFile >> EtaKPr; keysFile >> PhiKPr; keysFile >> DCAKPr; keysFile >> DCAxyKPr; keysFile >> DCAzKPr; keysFile >> PtKGl; keysFile >> PxKGl; keysFile >> PyKGl; keysFile >> PzKGl; keysFile >> EtaKGl; keysFile >> PhiKGl; keysFile >> DCAKGl; keysFile >> DCAxyKGl; keysFile >> DCAzKGl; keysFile >> GeantPi; keysFile >> StopRPi; keysFile >> KeyPi; keysFile >> RecoKeyPi; //cout << "Find Track RecoKeyPi-1 = " <> PtPiPr; keysFile >> PxPiPr; keysFile >> PyPiPr; keysFile >> PzPiPr; keysFile >> EtaPiPr; keysFile >> PhiPiPr; keysFile >> DCAPiPr; keysFile >> DCAxyPiPr; keysFile >> DCAzPiPr; keysFile >> PtPiGl; keysFile >> PxPiGl; keysFile >> PyPiGl; keysFile >> PzPiGl; keysFile >> EtaPiGl; keysFile >> PhiPiGl; keysFile >> DCAPiGl; keysFile >> DCAxyPiGl; keysFile >> DCAzPiGl; } cout << "FileId = " << FileId << " EventId[0] = " << EventId[0] << " mEventID = " << mEventId << endl; if(EventId[0] != mEventId || RunId[0] != mRunId) continue; foundEvent = 1; Double_t dEdxScale = 1; Int_t run = RunId[0]; #ifdef __OLD_DCA__ if (run < 10000) { dEdxScale = TMath::Exp(7.81779499999999961e-01); } else { if (run >= 5338005 && run <= 6177009) {// Run V, CuCu if (MuEvent_mRunInfo_mCenterOfMassEnergy[0] > 60 && MuEvent_mRunInfo_mCenterOfMassEnergy[0] < 70) {// 62 GeV // 76007 is cu-zdc-narrow (ZDC coincidence with 80 cm vertex cut), // 76011 is cu-bbc-narrow (BBC coincidence with vertex cut). // 76002 is cu-zdc-tacs (ZDC coincidence with no vertex cut.) if (run < 6069077 && ! (isTrigger(76002,mNTriggerId_mId) || isTrigger(76011,mNTriggerId_mId)) ) continue; if (run >= 6069077 && ! (isTrigger(76007,mNTriggerId_mId) || isTrigger(76011,mNTriggerId_mId)) ) continue; } else { /* 2005 CuCu http://www.star.bnl.gov/protected/common/common2005/trigger2005/CuFAQ.html Trigger Id Type Offline vertex cut Trigger setup 66007 Min Bias (cu-zdc-narrow) |vz|<30 cm cuProductionHighTower, cuProductionMinBias 66203 Barrel High Tower (cu-bemc-ht18)|vz|<30 cm cuProductionHighTower */ if (! (isTrigger(66007,mNTriggerId_mId) || isTrigger(66203,mNTriggerId_mId))) continue; } } } #endif /* 2007 AuAu http://www.star.bnl.gov/protected/common/common2007/trigger2007/triggers2007.html Production Trigger id's (remember that trigger id's < 1000 are test trigger id's with no guarantee as to sanity): Trigger Id Name First Run Last Run Description 200001 mb-vpd 8097121 8102062 ZDC coincidence+VPD cut at 5 cm. Note: the ZDC West was somewhat wide, may introduce inefficiencies. Killer bits on. 200003 mb-vpd 8103029 8113074 ZDC coincidence+VPD cut at 5 cm. Killer bits on. 200013 mb-vpd 8113077 8177059 ZDC coincidence+VPD cut at 5 cm. Killer bits off. 200020 mb30cm 8120052 8159044 ZDC coincidence+VPD cut at 30 cm. Only in 2007LowLuminosity trgsetupname. 200212 btag 8097121 8102029 ZDC coincidence + VPD cut at 5 cm + Barrel High Tower at threshold 1 (18, 4.1 GeV, accepting Et>4.3 GeV). Note: ZDC West wide, may lead to inefficiency. NOTE: this trigger id was reused for another meaning. 200213 btag 8103029 8113068 ZDC coincidence + VPD cut at 5 cm + Barrel High Tower at threshold 1 (18, 4.1 GeV, accepting Et>4.3 GeV). 200214 btag 8113102 8177038 ZDC coincidence + VPD cut at 5 cm + Barrel High Tower at threshold 1 (18, 4.1 GeV, accepting Et>4.3 GeV). ZDC killer bits off. 200601 L2-upsilon 8103029 8113068 L2 upsilon, based on VPD at 30 cm + L0 BHT1 (18, 4.08 GeV accepting Et>4.3 GeV) + L2 upsilon. No SSD required, but cannot happen at the same time as mb-vpd. 200602 L2-upsilon 8113102 8177038 L2 upsilon, based on ZDC + VPD at 30 cm + L0 BHT1 (18, 4.08 GeV accepting Et>4.3 GeV) + L2 upsilon. No SSD required, but cannot happen at the same time as mb-vpd. ZDC killer bits off 200610 upsilon-mb 8103029 8113068 Luminosity monitor for the L2-upsilon. ZDC coincidence + VPD at 30 cm. ZDC killer bits on. 200611 upsilon-mb 8113102 8177038 Luminosity monitor for the L2-upsilon. ZDC coincidence + VPD at 30 cm. ZDC killer bits off. 200620 L2-gamma 8109015 8113068 Tagger for express stream. Based on bht2-mb, with additional higher pt cluster cuts on cluster energy. 200621 L2-gamma 8113102 8177038 Tagger for express stream. Based on bht2-mb, with additional higher pt cluster cuts on cluster energy. */ #ifdef __TCFIT__ Double_t HZ = 0.000299792458 * mMagneticFieldZ[0]; #endif // if (NoTracks < 10) continue; //nPrim->Fill(NoPrimVertex); for (Int_t l = 0; l < NoPrimVertex; l++) { // if (! l) continue; //priVtxXY->Fill(PrimVertexX[l],PrimVertexY[l]); //priVtxZ->Fill(PrimVertexZ[l]); //priR->Fill(TMath::Sqrt(PrimVertexX[l]*PrimVertexX[l] +PrimVertexY[l]*PrimVertexY[l])); //priVtxSigmaZ->Fill(PrimVerSigZ[l]); //if (PrimVerSigZ[l] < 0 || PrimVerSigZ[l] > 0.080) continue; Double_t eta_max = - TMath::Log(TMath::Tan(0.5*TMath::ATan2(R, L-PrimVertexZ[l]))); Double_t eta_min = - TMath::Log(TMath::Tan(0.5*TMath::ATan2(R,-L-PrimVertexZ[l]))); Double_t vtx[3] = {PrimVertexX[l], PrimVertexY[l], PrimVertexZ[l]}; // pairs Int_t NoFSvtHits[2]; //[0] -> K, [1] -> pi Int_t NoFSsdHits[2]; Int_t NoFSvtSsdHits[2]; Int_t NoPSvtHits[2]; Int_t NoPSsdHits[2]; Int_t NoPSvtSsdHits[2]; Double_t dcaXY[2]; Double_t dcaZ[2]; Double_t sigmaXY[2]; Double_t sigmaZ[2]; TVector3 p[4]; TLorentzVector p4[4][2]; // K,pi; e,e; Int_t charge[2]; TVector3 dir[3]; TVector3 dcaG[2]; Double_t delPhi[4]; Double_t delTheta[4]; Double_t EffMass[4]; TLorentzVector PP[4]; Double_t cosL; Double_t cosP; Double_t sinP; #ifdef __TCFIT__ Double_t dL; Double_t eL; #endif Double_t slength; Double_t chisq; Double_t prob; Double_t dslength; //nTracks->Fill(NoTracks); Int_t NoTracksU = 0; Double_t sigmadEdx[2]; // //if (NoTracks > 50) continue; http://www.star.bnl.gov/protected/heavy/baumgart/D0_CuCu.html // Loop over primary tracks //k = RecoKeyK; //for (Int_t kg = 0; kg NoTracksGl) continue; //if (mNHitsFitGl[kg] < 15) continue; //if (pTGl[kg] < 0.1) continue; //if (pTGl[kg] < pTCut) continue; //if (EtaGl[kg] <= eta_min || EtaGl[kg] >= eta_max) continue; //if (dEdxTrackLength[kg] < 40.0) continue; sigmadEdx[0] = Bichsel::GetdEdxResolution(uTime[0], dEdxTrackLength[kg]); // if (mSigmaDcaDGl[kg] <= 0 || mSigmaDcaDGl[kg] > 1 || mSigmaDcaZGl[kg] <= 0 || mSigmaDcaZGl[kg] > 1) continue; NoFSvtHits[0] = (mNHitsFitInnerGl[kg] & 0x7); NoFSsdHits[0] = ((mNHitsFitInnerGl[kg] & 0x18) >> 3); NoFSvtSsdHits[0] = NoFSvtHits[0] + NoFSsdHits[0]; NoPSvtHits[0] = (mNHitsPossInnerGl[kg] & 0x7); NoPSsdHits[0] = ((mNHitsPossInnerGl[kg] & 0x18) >> 3); NoPSvtSsdHits[0] = NoPSvtHits[0] + NoPSsdHits[0]; p[0] = TVector3(pXGl[kg],pYGl[kg],pZGl[kg]); //<<< p[0] = TVector3(-5.5036, -3.6384,-1.2434); //p[2] = TVector3(pX[k],pY[k],pZ[k]); charge[0] = 0; if (QGl[kg] < 0) charge[0] = 1; dcaG[0] = TVector3(mDCAGlobal_mX1Gl[kg],mDCAGlobal_mX2Gl[kg],mDCAGlobal_mX3Gl[kg]); #ifdef __OLD_DCA__ dir[0] = p[0].Unit(); cosL = dir[0].Perp(); cosP = dir[0].x()/cosL; sinP = dir[0].y()/cosL; dcaXY[0] = sinP * dcaG[0].x() - cosP * dcaG[0].y(); // switched sign dcaZ[0] = dcaG[0].z()/(cosL*cosL); sigmaXY[0] = mSigmaDcaDGl[kg]; sigmaZ[0] = mSigmaDcaZGl[kg]/(cosL*cosL); // account bug in StMuTrack #else Int_t kgc = GlobalTracks_mIndex2Cov[kg]; static StDcaGeometry dcaGeometry; Float_t parsK[6] = { CovGlobTrack_mImp[kgc],CovGlobTrack_mZ[kgc],CovGlobTrack_mPsi[kgc], CovGlobTrack_mPti[kgc],CovGlobTrack_mTan[kgc],CovGlobTrack_mCurv[kgc]}; Float_t errsK[15] = { CovGlobTrack_mImpImp[kgc], CovGlobTrack_mZImp[kgc], CovGlobTrack_mZZ[kgc], CovGlobTrack_mPsiImp[kgc],CovGlobTrack_mPsiZ[kgc],CovGlobTrack_mPsiPsi[kgc], CovGlobTrack_mPtiImp[kgc],CovGlobTrack_mPtiZ[kgc],CovGlobTrack_mPtiPsi[kgc],CovGlobTrack_mPtiPti[kgc], CovGlobTrack_mTanImp[kgc],CovGlobTrack_mTanZ[kgc],CovGlobTrack_mTanPsi[kgc],CovGlobTrack_mTanPti[kgc],CovGlobTrack_mTanTan[kgc]}; dcaGeometry.set(parsK, errsK); THelixTrack thelix = dcaGeometry.thelix(); THelixTrack KaonTHelix = dcaGeometry.thelix();; StThreeVectorF momK = dcaGeometry.momentum(); StThreeVectorF ORIcopy = dcaGeometry.origin(); #ifdef __JB_HELICES__ //JB: for helices StThreeVectorF HP(hpX[kgc],hpY[kgc],hpZ[kgc]); StThreeVectorF ORI(OriginX[kgc],OriginY[kgc],OriginZ[kgc]); StPhysicalHelixD hel1(HP,ORI,hB[kgc]*kilogauss,QGl[kgc]); //JB #ifdef __JB_HELICES2__ StPhysicalHelixD hel1copy(momK,ORIcopy,hB[kgc]*kilogauss,QGl[kgc]); #endif #endif thelix.Move(thelix.Path(vtx)); dcaXY[0] = dcaGeometry.impact(); dcaZ[0] = dcaGeometry.z(); Double_t dcaEmx[3]; thelix.Dca(vtx,dcaXY[0],dcaZ[0],dcaEmx,2); sigmaXY[0] = 0.0; sigmaZ[0] = 0.0; if (dcaEmx[0] > 0) sigmaXY[0] = TMath::Sqrt(dcaEmx[0]); if (dcaEmx[2] > 0) sigmaZ[0] = TMath::Sqrt(dcaEmx[2]); #endif #if 1 sigmaXY[0] = TMath::Sqrt(sigmaXY[0]*sigmaXY[0] + PrimVerSigX[l]*PrimVerSigX[l] + PrimVerSigY[l]*PrimVerSigY[l]); sigmaZ[0] = TMath::Sqrt(sigmaZ[0]*sigmaZ[0] + PrimVerSigZ[l]*PrimVerSigZ[l]); #endif NoTracksU++; // Second Loop over primary tracks //i = RecoKeyPi; //for(Int_t ig = 0; ig NoTracksGl) continue; //if (pTGl[ig] < pTCut) continue; //if (EtaGl[ig] <= eta_min || EtaGl[ig] >= eta_max) continue; //if (mSigmaDcaDGl[ig] <= 0 || mSigmaDcaDGl[ig] > 1 || mSigmaDcaZGl[ig] <= 0 || mSigmaDcaZGl[ig] > 1) continue; //if (dEdxTrackLength[ig] < 40.0) continue; sigmadEdx[1] = Bichsel::GetdEdxResolution(uTime[0], dEdxTrackLength[ig]); NoFSvtHits[1] = (mNHitsFitInnerGl[ig] & 0x7); NoFSsdHits[1] = ((mNHitsFitInnerGl[ig] & 0x18) >> 3); NoFSvtSsdHits[1] = NoFSvtHits[1] + NoFSsdHits[1]; NoPSvtHits[1] = (mNHitsPossInnerGl[ig] & 0x7); NoPSsdHits[1] = ((mNHitsPossInnerGl[ig] & 0x18) >> 3); NoPSvtSsdHits[1] = NoPSvtHits[1] + NoPSsdHits[1]; p[1] = TVector3(pXGl[ig],pYGl[ig],pZGl[ig]); //p[3] = TVector3(pX[i],pY[i],pZ[i]); charge[1] = 0; if (QGl[ig] < 0) charge[1] = 1; dcaG[1] = TVector3(mDCAGlobal_mX1Gl[ig],mDCAGlobal_mX2Gl[ig],mDCAGlobal_mX3Gl[ig]); #ifdef __OLD_DCA__ dir[1] = p[1].Unit(); cosL = dir[1].Perp(); cosP = dir[1].x()/cosL; sinP = dir[1].y()/cosL; dcaXY[1] = sinP * dcaG[1].x() - cosP * dcaG[1].y(); // switched sign dcaZ[1] = dcaG[1].z()/(cosL*cosL); sigmaXY[1] = mSigmaDcaDGl[ig]; sigmaZ[1] = mSigmaDcaZGl[ig]/(cosL*cosL); // account bug in StMuTrack #else Int_t igc = GlobalTracks_mIndex2Cov[ig]; static StDcaGeometry dcaGeometryPi; Float_t parsPi[6] = { CovGlobTrack_mImp[igc],CovGlobTrack_mZ[igc],CovGlobTrack_mPsi[igc], CovGlobTrack_mPti[igc],CovGlobTrack_mTan[igc],CovGlobTrack_mCurv[igc]}; Float_t errsPi[15] = { CovGlobTrack_mImpImp[igc], CovGlobTrack_mZImp[igc], CovGlobTrack_mZZ[igc], CovGlobTrack_mPsiImp[igc],CovGlobTrack_mPsiZ[igc],CovGlobTrack_mPsiPsi[igc], CovGlobTrack_mPtiImp[igc],CovGlobTrack_mPtiZ[igc],CovGlobTrack_mPtiPsi[igc],CovGlobTrack_mPtiPti[igc], CovGlobTrack_mTanImp[igc],CovGlobTrack_mTanZ[igc],CovGlobTrack_mTanPsi[igc],CovGlobTrack_mTanPti[igc],CovGlobTrack_mTanTan[igc]}; dcaGeometryPi.set(parsPi, errsPi); THelixTrack thelixPi = dcaGeometryPi.thelix(); THelixTrack PionTHelix = dcaGeometryPi.thelix();; StThreeVectorF momPi = dcaGeometryPi.momentum(); StThreeVectorF ORIcopyPi = dcaGeometryPi.origin(); #ifdef __JB_HELICES__ //JB : for helices StThreeVectorF HP2(hpX[igc],hpY[igc],hpZ[igc]); StThreeVectorF ORI2(OriginX[igc],OriginY[igc],OriginZ[igc]); StPhysicalHelixD hel2(HP2,ORI2,hB[igc]*kilogauss,QGl[igc]); // JB #ifdef __JB_HELICES2__ StPhysicalHelixD hel2copy(momPi,ORIcopyPi,hB[igc]*kilogauss,QGl[igc]); #endif thelixPi.Move(thelixPi.Path(vtx)); //cout << "Moved Pi helix" << endl; dcaXY[1] = dcaGeometryPi.impact(); dcaZ[1] = dcaGeometryPi.z(); Double_t dcaEmx[3]; thelixPi.Dca(vtx,dcaXY[1],dcaZ[1],dcaEmx,2); sigmaXY[1] = 0; sigmaZ[1] = 0; if (dcaEmx[0] > 0) sigmaXY[1] = TMath::Sqrt(dcaEmx[0]); if (dcaEmx[2] > 0) sigmaZ[1] = TMath::Sqrt(dcaEmx[2]); #endif #endif #if 1 sigmaXY[1] = TMath::Sqrt(sigmaXY[1]*sigmaXY[1] + PrimVerSigX[l]*PrimVerSigX[l] + PrimVerSigY[l]*PrimVerSigY[l]); sigmaZ[1] = TMath::Sqrt(sigmaZ[1]*sigmaZ[1] + PrimVerSigZ[l]*PrimVerSigZ[l]); #endif TVector3 P(p[0]); P += p[1]; dir[2] = P.Unit(); /*if (k == 0 && i == 1 && EventId[0] == 4){ cout << "run = " << run << " k = " << k << " i = " << i << endl; cout << "pz0 = " << p[0].Z() << " pz1 = " << p[1].Z() << endl; cout << "phi0 = " << p[0].Phi() << " phi1 = " << p[1].Phi() << endl; cout << "pt0 = " << p[0].Perp() << " pt1 = " << p[1].Perp() << endl; cout << "eta0 = " << p[0].Eta() << " eta1 = " << p[1].Eta() << endl << endl; cout << "pz2 = " << P.Z() << endl; cout << "phi2 = " << P.Phi() << endl; cout << "pt2 = " << P.Perp() << endl; cout << "eta2 = " << P.Eta() << endl << endl; }*/ // Decay length linear approximation Double_t lTheta[3] = {p[0].Theta(), p[1].Theta(), P.Theta()}; Double_t lPhi[3] = {p[0].Phi(), p[1].Phi(), P.Phi()}; delPhi[0] = lPhi[0] - lPhi[2]; delPhi[1] = lPhi[1] - lPhi[2]; delTheta[0] = lTheta[0] - lTheta[2]; delTheta[1] = lTheta[1] - lTheta[2]; Double_t scaleXY = 1; Double_t scaleZ = 1; TRMatrix U(1,4, scaleXY*TMath::Sin(lTheta[2])*TMath::Sin(delPhi[0]), scaleZ*(TMath::Cos(lTheta[2]) - TMath::Sin(lTheta[2])/TMath::Tan(lTheta[0])*TMath::Cos(delPhi[0])), scaleXY*TMath::Sin(lTheta[2])*TMath::Sin(delPhi[1]), scaleZ*(TMath::Cos(lTheta[2]) - TMath::Sin(lTheta[2])/TMath::Tan(lTheta[1])*TMath::Cos(delPhi[1]))); TRSymMatrix W(4, 1./(sigmaXY[0]*sigmaXY[0]), 0. ,1./(sigmaZ[0]*sigmaZ[0]), 0. , 0.,1./(sigmaXY[1]*sigmaXY[1]), 0. , 0., 0., 1./(sigmaZ[1]*sigmaZ[1])); TRMatrix A(U, TRArray::kAxS, W); TRMatrix X(1,4, dcaXY[0], dcaZ[0], dcaXY[1], dcaZ[1]); TRMatrix B(A, TRArray::kAxBT, X); TRMatrix C(A, TRArray::kAxBT, U); //if (TMath::Abs(B[0]) < 1e-7) continue; slength = B[0]/C[0]; TRMatrix R(X); R = slength*U; TRSymMatrix RTWR(R,TRArray::kAxSxAT, W); chisq = RTWR[0]; prob = TMath::Prob(chisq,3); dslength = 1./TMath::Sqrt(C[0]); TVector3 XYZD(PrimVertexX[l],PrimVertexY[l],PrimVertexZ[l]); XYZD += slength*dir[2]; Float_t StopRD = 0.0; StopRD = XYZD.Perp(); //chisqR->Fill(chisq); //probR->Fill(prob); //if (slength > 3*dslength) { //decVtx->Fill(XYZD.X(),XYZD.Y()); //decR->Fill(TMath::Sqrt(XYZD.X()*XYZD.X() +XYZD.Y()*XYZD.Y())); //} #ifdef DEBUG cout << "U " << U << endl; cout << "W " << W << endl; cout << "A " << A << endl; cout << "B " << B << endl; cout << "C " << C << endl; cout << "X " << X << endl; cout << "R " << R << endl; cout << "RTWR " << RTWR << endl; cout << "chisq " << chisq << "\tProb = " << prob << endl; cout << "slength " << slength << " +/- " << dslength << " chisq " << chisq; if (dslength > 0) cout << " s.d.t " << slength/dslength; cout << endl; Double_t sav = 0, dsav = 0; for (Int_t i = 0; i < 4; i++) { Double_t uu = TMath::Abs(U[i]); Double_t ww = TMath::Sqrt(W(i,i)); Double_t uw2 = uu*uu*ww*ww; if (uu > 1e-7 && ww > 1e-7) { cout << "slength from [" << i << "] = " << X[i]/U[i] << " +/- " << 1./TMath::Sqrt(uw2) << endl; sav += X[i]/U[i]*uw2; dsav += uw2; } } if (dsav > 0) { sav /= dsav; dsav = 1./TMath::Sqrt(dsav); cout << " = " << sav << " +/- " << dsav << endl; } #endif /* DEBUG */ #ifdef __TCFIT__ //if (TMath::Abs(slength) > dLCut*dslength) { tc.Reset(); dat.Reset(); for (Int_t ip = 0; ip < 2; ip++) { dat.mTkBas[ip].Reset(); Double_t ptin = QGl[kg]/p[ip].Pt(); dat.mTkBas[ip].SetHz(1.); Double_t wk[9]; p[ip].GetXYZ(wk+3); dcaG[ip].GetXYZ(wk+0); THelixTrack th(wk+0,wk+3,ptin*HZ); th.Backward(); double s =th.Path(0.,0.); th.Move(s);th.Backward(); th.Eval(0.,wk,wk+3); TVector3 pos(wk); dir[ip] = TVector3(wk+3); dat.mTkBas[ip].Reset(); dat.mTkBas[ip].Set(pos,dir[ip],ptin); dat.mTkBas[ip].mass= (!ip)? amK:amPi; dat.mTkBas[ip].Update(); Float_t *errs = &errsK[0]; if (ip) errs = &errsPi[0]; for (Int_t n = 0; n < 5; n++) { for (Int_t m = 0; m <= n; m++) { Int_t nm = (n*(n+1))/2 + m; dat.mTEBas[ip].Set(n,m,errs[nm]); } } //dat.mTEBas[ip].Set(0,0,pow(sigmaXY[ip],2)); //dat.mTEBas[ip].Set(1,1,pow(sigmaZ[ip],2)); //Double_t dPhi = p[ip].Phi() - p[ip+2].Phi(); //Double_t dTanL = 1./TMath::Tan(p[ip].Theta()) - 1./TMath::Tan(p[ip+2].Theta()); //Double_t dpT = 1./p[ip].Pt() - 1./p[ip+2].Pt(); //dat.mTEBas[ip].Set(2,2,TMath::Max(1e-8,dPhi*dPhi)); // Phi //dat.mTEBas[ip].Set(3,3,TMath::Max(1e-8,dpT*dpT)); // 1./pT //dat.mTEBas[ip].Set(4,4,TMath::Max(1e-8,dTanL*dTanL)); // Tan dat.mTkBas[0].SetHz(mMagneticFieldZ[0]/4.98478); dat.mTkBas[1].SetHz(mMagneticFieldZ[0]/4.98478); } dat.Ready(); // dat.FixPar(TCFitV0::kPHI_0); // dat.FixPar(TCFitV0::kPHI_1); // dat.FixPar(TCFitV0::kPTIN_0); // dat.FixPar(TCFitV0::kPTIN_1); // dat.FixPar(TCFitV0::kTANL_0); // dat.FixPar(TCFitV0::kTANL_1); dat.FixPar(TCFitV0::kCNRJ); tc.SetMaxIter(10); //cout << "Made it Here! " << endl; if (tc.Fit()) { #ifdef DEBUG cout << "tc.Fit fails" << endl; #endif continue; } dL = dat.GetPar(TCFitV0::kLEN_2); eL = sqrt(dat.ErMx(TCFitV0::kLEN_2,TCFitV0::kLEN_2)); chisq = dat.GetFcn(); prob = TMath::Prob(chisq,dat.GetNDF()); #ifdef DEBUG cout << "dL = " << dL << " +/- " << eL << " chisq = " << chisq << " NDF = " << dat.GetNDF() << "\tprob = " << prob << endl; #endif //SvsL->Fill(dL,slength); slength = dL; dslength = eL; TLorentzVector ParentMomTCFit = dat.mTkFit[0].P4()+dat.mTkFit[1].P4(); Double_t decayTCFitX = dL*ParentMomTCFit.Px()/ParentMomTCFit.P(); Double_t decayTCFitY = dL*ParentMomTCFit.Py()/ParentMomTCFit.P(); Double_t decayTCFitZ = dL*ParentMomTCFit.Pz()/ParentMomTCFit.P(); TVector3 decayTCFit; decayTCFit.SetXYZ(decayTCFitX,decayTCFitY,decayTCFitZ); Double_t decayLengthTCFit = decayTCFit.Mag(); Double_t stopRDTCFit = decayTCFit.Perp(); //} #endif //if (TMath::Abs(slength) > 0.20) continue; // Int_t ki[2][2] = {{kg,ig},{ig,kg}}; // Int_t kip[2][2] = {{0,1},{1,0}}; Int_t c = charge[0] + 2*charge[1]; Double_t dEdx[2] = {dEdxScale*GdEdx[kg],dEdxScale*GdEdx[ig]}; Double_t am[4] = {amK, amPi, amP, ame}; Double_t pId[2][4]; for (Int_t t = 0; t < 2; t++) { for (Int_t h = 0; h < 4; h++) { Double_t bg = p[t].Mag()/am[h]; Double_t bg10 = TMath::Log10(bg); Double_t dEdxP = 1.e-6*m_Bichsel->GetI70M(bg10,1,0); pId[t][h] = TMath::Log(dEdx[t]/dEdxP)/sigmadEdx[t]; } } #if 0 Double_t pIdO[2][4] = { {NSigmaKaonGl[kg]/__SIGMA_SCALE__, NSigmaPionGl[kg]/__SIGMA_SCALE__, NSigmaProtonGl[kg]/__SIGMA_SCALE__, NSigmaElectronGl[kg]/__SIGMA_SCALE__}, {NSigmaKaonGl[ig]/__SIGMA_SCALE__, NSigmaPionGl[ig]/__SIGMA_SCALE__, NSigmaProtonGl[ig]/__SIGMA_SCALE__, NSigmaElectronGl[ig]/__SIGMA_SCALE__} }; #endif Bool_t KaonPiD = TMath::Abs(pId[0][0]) < 2.0 && TMath::Abs(pId[0][1]) > 2.0 && TMath::Abs(pId[0][2]) > 2.0;// && TMath::Abs(pId[0][3]) > 2.0; Bool_t PionPiD = TMath::Abs(pId[1][1]) < 2.0 && TMath::Abs(pId[1][0]) > 2.0 && TMath::Abs(pId[1][2]) > 2.0;// && TMath::Abs(pId[1][3]) > 2.0; //pIdKpi[c]->Fill(pId[1][1],pId[0][0]); //TdEdx->Fill(TMath::Log10(p[0].Mag()), TMath::Log10(dEdx[0])+6.); //TdEdx->Fill(TMath::Log10(p[1].Mag()), TMath::Log10(dEdx[1])+6.); Bool_t KpiPiD = KaonPiD && PionPiD; //if (KpiPiD) { //pIdKpiC[c]->Fill(pId[1][1],pId[0][0]); //TdEdxKaonC->Fill(TMath::Log10(p[0].Mag()), TMath::Log10(dEdx[0])+6.); //TdEdxPionC->Fill(TMath::Log10(p[1].Mag()), TMath::Log10(dEdx[1])+6.); //} p4[0][0].SetVectMag(p[0],amK); p4[1][0].SetVectMag(p[1],amPi); //p4[0][1].SetVectMag(p[2],amK); //p4[1][1].SetVectMag(p[3],amPi); PP[0] = p4[0][0]; PP[0] += p4[1][0]; //} /* The Gottfried - Jackson (GJ) frame is a rest frame of the Kpi system in which the z-axis is in the direction of the beam momentum and the y-axis is in the direction of the vector cross-product of the target and recoil momenta. Gottfried - Jackson: cosTheta is projection in GJ frame momentum of K on direction of Kpi in original frame */ TVector3 bF = PP[0].BoostVector(); TVector3 b(-bF.X(),-bF.Y(),-bF.Z()); TLorentzVector Kl(p4[0][0]); Kl.Boost(b); cout << "Klx = " << Kl.X() << endl; TVector3 dother(Kl.Vect()); TVector3 mother(PP[0].Vect()); Double_t cosTheta_GJ_K = dother.Dot(mother)/(dother.Mag()*mother.Mag()); EffMass[0] = PP[0].M(); //PP[1] = p4[0][1]; //PP[1] += p4[1][1]; //EffMass[1] = PP[1].M(); for (Int_t f = 0; f < NF; f++) { TLorentzVector Pfake(p4[0][0]); TLorentzVector pfake(p4[1][0]); //Double_t rot = TMath::Pi()/180*(150 + 60/(NF - 1)*f); // 150 <= 210 with 5 degree step Double_t rot = TMath::Pi(); // Rotate by 180 degrees pfake.RotateZ(rot); Pfake += pfake; if (2*f + 1 == NF) PP[2] = Pfake; EffMass[2] = Pfake.M(); } #ifdef DEBUG cout << "MKpi_GL\t" << EffMass[0] << "\tMKpi_Pr\t" << EffMass[1] << "\tMKpi_Fake\t" << EffMass[2] << "\tCos(Theta_GJ)\t" << cosTheta_GJ_K << endl;; #endif /* DEBUG */ #ifdef __JB_HELICES__ // use THelixTrack #ifdef __JB_HELICES2__ double d1=0,d2=0,x1[3],x2[3]; for(int i=0;i<3;i++) { x1[i] = 0.0; x2[i] = 0.0; } d1 = KaonTHelix.Path(PionTHelix,&d2); KaonTHelix.Eval(d1,x1); PionTHelix.Eval(d2,x2); Float_t deltaX = x1[0] - x2[0]; Float_t deltaY = x1[1] - x2[1]; Float_t deltaZ = x1[2] - x2[2]; StThreeVectorD X1T(x1[0],x1[1],x1[2]); StThreeVectorD X2T(x2[0],x2[1],x2[2]); StThreeVectorD MidPointT = (X1T+X2T)/2; StThreeVectorD mPVT(PrimVertexX[l],PrimVertexY[l],PrimVertexZ[l]); StThreeVectorD decayT = MidPointT - mPVT ; #endif //JB : for helices pairD pathLengthtempo = hel1.pathLengths(hel2); StThreeVectorD hPos = hel1.at(pathLengthtempo.first); StThreeVectorD ePos = hel2.at(pathLengthtempo.second); StThreeVectorD vecDca = ePos-hPos; //float trackDca = vecDca.mag(); //distance between the 2 tracks StThreeVectorD X1 = hel1.at(pathLengthtempo.first); StThreeVectorD X2 = hel2.at(pathLengthtempo.second); StThreeVectorD MidPoint = (X1+X2)/2; StThreeVectorD mPV(PrimVertexX[l],PrimVertexY[l],PrimVertexZ[l]); StThreeVectorD decay = MidPoint - mPV ; StThreeVectorD P0(pXGl[kg],pYGl[kg],pZGl[kg]); // first loop StThreeVectorD P1(pXGl[ig],pYGl[ig],pZGl[ig]); // second loop #ifdef __JB_HELICES2__ // test helix at (0,0,0) pairD pathLengthcopy = hel1copy.pathLengths(hel2copy); StThreeVectorD hPoscopy = hel1copy.at(pathLengthcopy.first); StThreeVectorD ePoscopy = hel2copy.at(pathLengthcopy.second); StThreeVectorD vecDcacopy = ePoscopy-hPoscopy; StThreeVectorD X1C = hel1copy.at(pathLengthcopy.first); StThreeVectorD X2C = hel2copy.at(pathLengthcopy.second); StThreeVectorD MidPointC = (X1C+X2C)/2; StThreeVectorD decayC = MidPointC - mPV; #endif StThreeVectorD PTOT = P0+P1; //sum of momentum of daughters StThreeVectorD CROSS = decay.cross(PTOT); Double_t DCAXY_p = TMath::Sqrt(CROSS.x()*CROSS.x() + CROSS.y()*CROSS.y()); DCAXY_p = DCAXY_p/PTOT.mag(); Double_t DCAZ_p = CROSS.z()/PTOT.mag(); //JB : for helices #endif //if(EventId[0] == 81 && NumFiles ==0) cout << "event 81 found D0" << endl; //if (KpiPiD){ //D0found++; //cout << "D0 found NumD0's = " << D0found << " in NevProc = " << NevProc << endl; //cout << "EventID = " << EventId[0] << " RunId = " << RunId[0] << " k = " << k << " i = " << i << endl; //cout << "kg = " << kg << " ig = " << ig << " pt1 = " << p4[0][0].Pt() << " eta1 = " << p4[0][0].Eta() << " pz1 = " << p4[0][0].Pz() << endl; //cout << "pt2 = " << p4[1][0].Pt() << " eta2 = " << p4[1][0].Eta() << " pz2 = " << p4[1][0].Pz() << endl; //cout << "ptP = " << PP[0].Pt() << " etaP = " << PP[0].Eta() << " pzP = " << PP[0].Pz() << " mP = " << PP[0].M() << endl << endl << endl; cout << "FileId = " << fileIdIntMuDst << endl; //cout << "VtxDx = " << XYZD.X() << " VtxDy = " << XYZD.Y() << " VtxDz = " << XYZD.Z() << endl; cout << "DCAxyK = " << dcaXY[0] << " DCAzK = " << dcaZ[0] << endl; cout << "DCAxyPi = " << dcaXY[1] << " DCAzPi = " << dcaZ[1] << endl; cout << "decayTCFitX = " << decayTCFit.X() << " decayTCFitY = " << decayTCFit.Y() << " decayTCFitZ = " << decayTCFit.Z() << endl; cout << "decayHelixX = " << decay.x() << " decayHelixY = " << decay.y() << " decayHelixZ = " << decay.z() << endl; D0found++; //if (PP[0].M() > 4){ //cout << "FileId = " << NumFiles << " EventId = " << EventId[0] << " M = " << PP[0].M() << endl << endl; //cout << "EventID = " << EventId[0] << " RunId = " << RunId[0] << " k = " << k << " i = " << i << endl; //cout << "kg = " << kg << " ig = " << ig << " pt1 = " << p4[0][0].Pt() << " eta1 = " << p4[0][0].Eta() << " pz1 = " << p4[0][0].Pz() << endl; //cout << "ptKGl = " << PtKGl << endl; //cout << "pt2 = " << p4[1][0].Pt() << " eta2 = " << p4[1][0].Eta() << " pz2 = " << p4[1][0].Pz() << endl; //cout << "ptPiGl = " << PtPiGl << endl; //cout << "ptP = " << PP[0].Pt() << " etaP = " << PP[0].Eta() << " pzP = " << PP[0].Pz() << " mP = " << PP[0].M() << endl << endl << endl; //} Int_t NoTracks = 0; //To Deal with lack of primary tracks TVector3 decayOrig(XYZD.X()-PrimVertexX[l],XYZD.Y()-PrimVertexY[l],XYZD.Z()-PrimVertexZ[l]); Float_t decayLength = TMath::Sqrt(decayOrig.X()*decayOrig.X()+decayOrig.Y()*decayOrig.Y()+decayOrig.Z()*decayOrig.Z()); Float_t decayLengthHel1 = TMath::Sqrt(decay.x()*decay.x()+decay.y()*decay.y()+decay.z()*decay.z()); Float_t stopRDHel1 = TMath::Sqrt(decay.x()*decay.x()+decay.y()*decay.y()); Float_t decayLengthHel2 = TMath::Sqrt(decayC.x()*decayC.x()+decayC.y()*decayC.y()+decayC.z()*decayC.z()); Float_t stopRDHel2 = TMath::Sqrt(decayC.x()*decayC.x()+decayC.y()*decayC.y()); fprintf(muDstData,"%i \t %i \t %i \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %i \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %i \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %i \t %i \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \t %.4f \n", fileIdIntMuDst,EventId[0],RunId[0],PrimVertexX[l],PrimVertexY[l],PrimVertexZ[l], PP[0].Pt(),PP[0].Px(),PP[0].Py(),PP[0].Pz(),PP[0].Eta(),PP[0].M(),XYZD.X(),XYZD.Y(),XYZD.Z(),cosTheta_GJ_K, kg+1,p4[0][0].Pt(),p4[0][0].Pz(),p4[0][0].Eta(),dcaXY[0],dcaZ[0], ig+1,p4[1][0].Pt(),p4[1][0].Pz(),p4[1][0].Eta(),dcaXY[1],dcaZ[1],NoTracksGl,NoTracks,decayLength,sigmaXY[0],sigmaXY[1],slength,p4[0][0].Px(),p4[0][0].Py(),p4[1][0].Px(),p4[1][0].Py(), StopRD,decayOrig.X(),decayOrig.Y(),decayOrig.Z(),decay.x(),decay.y(),decay.z(),decayLengthHel1,stopRDHel1,decayC.x(),decayC.y(),decayC.z(),decayLengthHel2,stopRDHel2,decayTCFit.X(),decayTCFit.Y(),decayTCFit.Z(),decayLengthTCFit,stopRDTCFit); tuple[0] = (Float_t)fileIdIntMuDst; tuple[1] = (Float_t)EventId[0]; tuple[2] = (Float_t)RunId[0]; tuple[3] = PrimVertexX[l]; tuple[4] = PrimVertexY[l]; tuple[5] = PrimVertexZ[l]; tuple[6] = PP[0].Pt(); tuple[7] = PP[0].Px(); tuple[8] = PP[0].Py(); tuple[9] = PP[0].Pz(); tuple[10] = PP[0].Eta(); tuple[11] = PP[0].M(); tuple[12] = XYZD.X(); tuple[13] = XYZD.Y(); tuple[14] = XYZD.Z(); tuple[15] = cosTheta_GJ_K; tuple[16] = (Float_t) kg; tuple[17] = p4[0][0].Pt(); tuple[18] = p4[0][0].Pz(); tuple[19] = p4[0][0].Eta(); tuple[20] = dcaXY[0]; tuple[21] = dcaZ[0]; tuple[22] = sigmaXY[0]; tuple[23] = sigmaZ[0]; tuple[24] = (Float_t) ig; tuple[25] = p4[1][0].Pt(); tuple[26] = p4[1][0].Pz(); tuple[27] = p4[1][0].Eta(); tuple[28] = dcaXY[1]; tuple[29] = dcaZ[1]; tuple[30] = sigmaXY[1]; tuple[31] = sigmaZ[1]; tuple[32] = (Float_t)NoTracksGl; tuple[33] = (Float_t)NoTracks; //tuple[34] = sigmaXY[0]; //tuple[35] = sigmaXY[1]; tuple[34] = PP[0].Px()-PxDMc; //cout << "PxDMuKpi = " << PP[0].Px() << " PxDMc = " << PxDMc << endl; //if(TMath::Abs(PP[0].Px()-PxDMc) > 0.5) cout << "PxDMuKpi - PxDMc = " << PP[0].Px()-PxDMc << endl; tuple[35] = PP[0].Py()-PyDMc; tuple[36] = PP[0].Pz()-PzDMc; tuple[37] = PP[0].Eta()-EtaDMc; tuple[38] = PP[0].Phi()-PhiDMc; tuple[39] = decayLength; tuple[40] = slength; //tuple[40] = decayLengthD - decayLength; nTuple->Fill(tuple); //if ( mKpiMin < mKpiMax) { //if (EffMass[0] < mKpiMin || EffMass[0] > mKpiMax) continue; //} //pTKpi[c]->Fill(p[1].Perp(),p[0].Perp()); //pKpi[c]->Fill(p[1].Mag(),p[0].Mag()); //if (p[0].Mag() < 0.5 || p[1].Mag() < 0.5) continue; // || p[0].Mag()*p[1].Mag() < 0.65) continue; //pTKpiC[c]->Fill(p[1].Perp(),p[0].Perp()); //pKpiC[c]->Fill(p[1].Mag(),p[0].Mag()); //CosTGJKpi[c]->Fill(cosTheta_GJ_K); //DCAxyKpi[c]->Fill(dcaXY[1],dcaXY[0]); //DCAzKpi[c]->Fill(dcaZ[1],dcaZ[0]); //dEdxKpi[c]->Fill(dEdx[1],dEdx[0]); /*for (Int_t s = 0; s < Nsys; s++) { Double_t pT = PP[s].Perp(); Int_t nz = 1; if (TMath::Abs(PrimVertexZ[l]) < 10 ) nz = NZ; else if (TMath::Abs(PrimVertexZ[l]) < 20 ) nz = NZ - 1; else if (TMath::Abs(PrimVertexZ[l]) < 30 ) nz = NZ - 2; Int_t neta = -1, npT = -1; for (npT = NpT-1; npT > 0; npT--) { //cout << npT << "\tpT " << pT << " pTmin " << pTmin[npT] << endl; if (pT > pTmin[npT]) break; } if (npT >= NpT) npT = 0; // for (npT = 1; npT < NpT; npT++) {if (pT < pTmin[npT]) break;} Double_t eta = TMath::Abs(PP[s].Eta()); // for (neta = 1; neta < Neta; neta++) {if (eta < etamin[neta]) break;} for (neta = 1; neta < Neta; neta++) {if (eta < etamin[neta]) break;} if (neta >= Neta) neta = 0; Double_t STD = -99; if (dslength > 0) STD = slength/dslength; #if 0 // 0 DL,1 SL, 2 EM, 3 RD, 4 C2, 5 DS Double_t Vars[6] = {TMath::Abs(slength), TMath::Abs(STD), EffMass[s], -1, chisq, dslength}; if (EffMass[1] < 0.080) Vars[5] = XYZD.Perp(); #else // 0 DL,1 SL, 2 EM, EMVx , EMF Double_t Vars[5] = {TMath::Abs(slength), STD, EffMass[0], EffMass[1], EffMass[2]}; #endif // Bool_t LDCAcut = prob > 1e-2 && (dslength > 0 && TMath::Abs(slength) < DcaCut && // TMath::Abs(dcaXY[0]) < DcaCut && TMath::Abs(dcaXY[1]) < DcaCut && // TMath::Abs(dcaZ[0]) < DcaCut && TMath::Abs(dcaZ[1]) < DcaCut); /* for (Int_t z = 0; z < nz; z++) { for (Int_t lk = 0; lk < NL; lk++) { if (lk == 0 || // ((prob > 1e-2 && dslength > 0) && (LDCAcut && (lk == 1 || lk == 2 && slength > 0 || lk == 3 && slength < 0 || lk == 4 && STD > 1 || lk == 5 && STD < -1 || lk == 6 && STD > 2 || lk == 7 && STD < -2 || lk == 8 && STD > 3 || lk == 9 && STD < -3))) { for (Int_t mGJ = 0; mGJ < NGJ; mGJ++) { if (mGJ == 0 || mGJ == 1 && TMath::Abs(cosTheta_GJ_K) <= 0.6 || mGJ == 2 && TMath::Abs(cosTheta_GJ_K) > 0.6) { for (Int_t cut = 0; cut < Ncut; cut++) { if (cut == 0 || cut == 1 && NoFSvtSsdHits[0] > 0 && NoFSvtSsdHits[1] > 0 || cut == 4 && NoFSvtSsdHits[0] > 1 && NoFSvtSsdHits[1] > 1 || (cut == 2 || cut == 3 && NoFSvtSsdHits[0] > 0 && NoFSvtSsdHits[1] > 0 || cut == 5 && NoFSvtSsdHits[0] > 1 && NoFSvtSsdHits[1] > 1) && KpiPiD) { for (Int_t kpT = 0; kpT <= npT; kpT++) { for (Int_t keta = 0; keta <= neta; keta++) { for (Int_t t = 0; t < NT-1; t++) { //hists[s][c][t][z][lk][mGJ][cut][kpT][keta]->Fill(Vars[t]); for (Int_t f = 0; f < NF; f++) { //hists[s][c][NT-1][z][lk][mGJ][cut][kpT][keta]->Fill(EffMass[2+f]); } } } } } } } } } } }*/ //} } } } } //nTracksU->Fill(NoTracksU); } if (NevProc%200 == 1) cout << NevProc << "\tevents processed so far" << endl; //if (NevProc > 5) break; if(foundEvent){ keysFile >> FileId; keysFile >> mEventId; //cout << "Find mEventId = " << mEventId << endl; keysFile >> mRunId; keysFile >> mVertexX; keysFile >> mVertexY; keysFile >> mVertexZ; keysFile >> StopRDMc; keysFile >> PtDMc; keysFile >> PxDMc; keysFile >> PyDMc; keysFile >> PzDMc; keysFile >> EtaDMc; keysFile >> PhiDMc; keysFile >> GeantK; keysFile >> StopRK; keysFile >> KeyK; keysFile >> RecoKeyK; //cout << "Find Track RecoKeyK-1 = " << RecoKeyK-1 << endl; keysFile >> PtKPr; keysFile >> PxKPr; keysFile >> PyKPr; keysFile >> PzKPr; keysFile >> EtaKPr; keysFile >> PhiKPr; keysFile >> DCAKPr; keysFile >> DCAxyKPr; keysFile >> DCAzKPr; keysFile >> PtKGl; keysFile >> PxKGl; keysFile >> PyKGl; keysFile >> PzKGl; keysFile >> EtaKGl; keysFile >> PhiKGl; keysFile >> DCAKGl; keysFile >> DCAxyKGl; keysFile >> DCAzKGl; keysFile >> GeantPi; keysFile >> StopRPi; keysFile >> KeyPi; keysFile >> RecoKeyPi; //cout << "Find Track RecoKeyPi-1 = " <> PtPiPr; keysFile >> PxPiPr; keysFile >> PyPiPr; keysFile >> PzPiPr; keysFile >> EtaPiPr; keysFile >> PhiPiPr; keysFile >> DCAPiPr; keysFile >> DCAxyPiPr; keysFile >> DCAzPiPr; keysFile >> PtPiGl; keysFile >> PxPiGl; keysFile >> PyPiGl; keysFile >> PzPiGl; keysFile >> EtaPiGl; keysFile >> PhiPiGl; keysFile >> DCAPiGl; keysFile >> DCAxyPiGl; keysFile >> DCAzPiGl; foundEvent = 0; } //if(mEventId == EventId[0] && EventId[0] < nEventsMax){ // nPairs++; // foundEvent = 1; //} //} } cout << "***** I have found " << D0found << " D0s. *****" << endl; } inFile.close(); keysFile.close(); } fclose(muDstData); fOut->Write(); delete fOut; } /* EMKpiNP 0.108 EMKpiNPLp 0.116 EMKpiNPLn 0.232 EMKpiNPLp1 0.134 EMKpiNPLn1 0.748 EMKpiNPLp2 0.157 EMKpiNPLp2 1.711 EMKpiNPLp3 0.381 EMKpiNPLn3 0.772 EMKpiNPpT05 0.081 EMKpiNPLppT05 0.090 EMKpiNPpT10 0.119 EMKpiNPLppT10 0.110 EMKpiNPpT20 0.096 EMKpiNPLppT20 0.099 EMKpiNPLnpT20 0.217 EMKpiNPGJPpT20 0.172 EMKpiNPLpGJPpT20 0.136 EMKpiNPLnGJPpT20 0.128 EMKpiNPGJNpT20 0.123 EMKpiNPLpGJNpT20 0.143 EMKpiNPLnGJNpT20 0.180 EMKpiNPLp1pT20 0.104 EMKpiNPLn1pT20 0.515 EMKpiNPLp2GJPpT20 0.144 EMKpiNPLp2GJNpT20 0.254 EMKpiNPLp2pT20 0.130 EMKpiNPLn2pT20 1.173 EMKpiNPLp3pT20 0.507 EMKpiNPLp3pT20 9.999 01/20/08 KpiD0Mix #std cucu200c 63.4 M, |z| < 10 -> 20M) EMKpiNP 12.3742 EMKpiNPLp 10.622 EMKpiNPLp1 8.28414 EMKpiNPLp2 6.77897 EMKpiNPLp3 2.39407 EMKpiNPLn 7.284 EMKpiNPLn1 4.77356 EMKpiNPLn2 3.43843 EMKpiNPLn3 1.57657 EMKpiNPdEdx 12.4304 EMKpiNPGJPdEdx 7.80157 EMKpiNPGJNdEdx 10.013 EMKpiNP 12.3742 EMKpiNPpT05 10.9887 EMKpiNPpT10 8.85068 EMKpiNPpT20 10.2894 EMKpiNPdEdxpT20 8.66682 EMKpiNPGJPdEdxpT20 5.69601 EMKpiNPGJNdEdxpT20 6.42735 EMKpiNPLpGJNdEdxpT20 3.14536 EMKpiNPLnGJNdEdxpT20 5.04419 EMKpiNPpT50 3.461 EMKpiPN No signal EMKpiNPpT10 4.01422 EMKpiPNpT10 3.42876 */