* provided "as is" without express or implied warranty. *\r
**************************************************************************/\r
\r
-////////////////////////////////////////////////////////////////////////////////\r
+////////////////////////////////////////////////////////////////////////////\r
//\r
-// This class is used to perform femtoscopic analysis on K0s particles, which\r
-// are reconstructed using the AliAODv0 class. \r
+// This class is used to perform femtoscopic analysis on K0s particles, \r
+// which are reconstructed using the AliAODv0 class. \r
//\r
// authors: Matthew Steinpreis (matthew.steinpreis@cern.ch)\r
// \r
-// Major changes:\r
+// Change log:\r
// - TOF mismatch function calls changed (4/18/13)\r
-// - added minimum decay length cut (3/28/13)\r
-//\r
-// Minor changes:\r
+// - added minimum decay length cut (rarely used though) (3/28/13)\r
// - K0 multiplicity histogram now filled with "unskippedCount" instead\r
-// of k0Count (which included skipped k0s with shared daughters) (3/25/13)\r
+// of k0Count (which included skipped k0s with shared daughters) \r
+// (3/25/13)\r
// - added hists for 3D mom. in LF and PRF (3/28/13) \r
// - changed calling of PIDResponse (should be same actions) (3/28/13)\r
// - keep "side" K0s for mass plot (4/18)\r
// - tweaked loading and skipping appropriately\r
// - use merit test to skip sides (against good and side K0s)\r
// - a good K0 cant be skipped by a side\r
+// - moved TPC propagation (via Hans' method) up to v0 level, which now\r
+// uses an AliAODTrack(AliVTrack) instead of AliESDtrack (5/31/13)\r
+// - added primary vertex subtraction in TPC propagation (5/31/13)\r
+// - removed all instances of AliESDtrack usage (5/31/13)\r
+// - removed old separation method/histograms (5/31/13)\r
+// - tidied up LCMS boost (6/10/13)\r
+// - added new boosting prescription, get q out-side-long for LCMS and PRF (6/24/13)\r
+// - added histograms and values for LCMS momenta (for simulation)\r
+// - added random particle order switch in correlations (9/09/13)\r
+// - added more bins for 3D OSL analysis (9/19/13)\r
+// - added merit cut choice, pass as argument (10/16/13)\r
+// - 1-mass, 2-v0dca, 3-dddca, 4-combination (used to be v0dca)\r
+// - added passable argument for two-track minimum separation (10/16/13)\r
////////////////////////////////////////////////////////////////////////////////\r
\r
\r
-\r
+ \r
#include <iostream>\r
#include <math.h>\r
#include "TMath.h"\r
fOnlineCase(kTRUE),\r
fMeritCase(kTRUE),\r
fMinDecayLength(0.0),\r
+ fMeritCutChoice(0),\r
+ fMinSep(0.0),\r
fEventCount(0),\r
fEC(0x0),\r
fEvt(0X0),\r
fName(0x0),\r
fAOD(0x0),\r
fOutputList(0x0),\r
- fPidAOD(0x0),\r
- fPosDaughter1(0x0), \r
- fPosDaughter2(0x0),\r
- fNegDaughter1(0x0),\r
- fNegDaughter2(0x0)\r
+ fPidAOD(0x0)\r
{\r
}\r
//________________________________________________________________________\r
-AliFemtoK0Analysis::AliFemtoK0Analysis(const char *name, bool FieldPositive, bool OnlineCase, bool MeritCase, float MinDL) \r
+AliFemtoK0Analysis::AliFemtoK0Analysis(const char *name, bool FieldPositive, bool OnlineCase, bool MeritCase, float MinDL, int MeritCutChoice, float MinSep) \r
: AliAnalysisTaskSE(name), \r
fFieldPos(FieldPositive),\r
fOnlineCase(OnlineCase),\r
fMeritCase(MeritCase),\r
fMinDecayLength(MinDL),\r
+ fMeritCutChoice(MeritCutChoice),\r
+ fMinSep(MinSep),\r
fEventCount(0),\r
fEC(0x0),\r
fEvt(0X0),\r
fName(name),\r
fAOD(0x0),\r
fOutputList(0x0),\r
- fPidAOD(0x0),\r
- fPosDaughter1(0x0), \r
- fPosDaughter2(0x0),\r
- fNegDaughter1(0x0),\r
- fNegDaughter2(0x0)\r
+ fPidAOD(0x0)\r
{\r
//main constructor\r
fFieldPos = FieldPositive;\r
fOnlineCase = OnlineCase;\r
fMeritCase = MeritCase;\r
fMinDecayLength = MinDL;\r
+ fMeritCutChoice = MeritCutChoice;\r
+ fMinSep = MinSep;\r
\r
// Define output slots here \r
// Output slot #1\r
fOnlineCase(obj.fOnlineCase),\r
fMeritCase(obj.fMeritCase),\r
fMinDecayLength(obj.fMinDecayLength),\r
+ fMeritCutChoice(obj.fMeritCutChoice),\r
+ fMinSep(obj.fMinSep),\r
fEventCount(obj.fEventCount),\r
fEC(obj.fEC),\r
fEvt(obj.fEvt),\r
fName(obj.fName),\r
fAOD(obj.fAOD),\r
fOutputList(obj.fOutputList),\r
- fPidAOD(obj.fPidAOD),\r
- fPosDaughter1(obj.fPosDaughter1), \r
- fPosDaughter2(obj.fPosDaughter2),\r
- fNegDaughter1(obj.fNegDaughter1),\r
- fNegDaughter2(obj.fNegDaughter2)\r
+ fPidAOD(obj.fPidAOD)\r
{\r
}\r
//________________________________________________________________________\r
fOnlineCase = obj.fOnlineCase;\r
fMeritCase = obj.fMeritCase;\r
fMinDecayLength= obj.fMinDecayLength;\r
+ fMeritCutChoice= obj.fMeritCutChoice;\r
+ fMinSep = obj.fMinSep;\r
fEventCount = obj.fEventCount;\r
fEC = obj.fEC;\r
fEvt = obj.fEvt;\r
fAOD = obj.fAOD;\r
fOutputList = obj.fOutputList;\r
fPidAOD = obj.fPidAOD;\r
- fPosDaughter1 = obj.fPosDaughter1; \r
- fPosDaughter2 = obj.fPosDaughter2;\r
- fNegDaughter1 = obj.fNegDaughter1;\r
- fNegDaughter2 = obj.fNegDaughter2;\r
\r
return (*this);\r
}\r
if(fAOD) delete fAOD;\r
if(fOutputList) delete fOutputList;\r
if(fPidAOD) delete fPidAOD;\r
- if(fPosDaughter1) delete fPosDaughter1;\r
- if(fPosDaughter2) delete fPosDaughter2;\r
- if(fNegDaughter1) delete fNegDaughter1;\r
- if(fNegDaughter2) delete fNegDaughter2;\r
}\r
//________________________________________________________________________\r
void AliFemtoK0Analysis::MyInit()\r
AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
fPidAOD = aodH->GetPIDResponse();\r
\r
- fPosDaughter1 = new AliESDtrack();\r
- fPosDaughter2 = new AliESDtrack();\r
- fNegDaughter1 = new AliESDtrack();\r
- fNegDaughter2 = new AliESDtrack();\r
-\r
fRandomNumber = new TRandom3(); //for 3D, random sign switching\r
fRandomNumber->SetSeed(0);\r
\r
TH2F* fHistKtK0 = new TH2F("fHistKtK0", "Kt distribution of K0 pairs", kCentBins, .5, kCentBins+.5, 300, 0., 3.);\r
fOutputList->Add(fHistKtK0);\r
\r
+ TH1F* fHistPx = new TH1F("fHistPx","",200,0,2);\r
+ TH1F* fHistPy = new TH1F("fHistPy","",200,0,2);\r
+ TH1F* fHistPz = new TH1F("fHistPz","",200,0,2);\r
+ TH1F* fHistPxCM = new TH1F("fHistPxCM","",200,0,2);\r
+ TH1F* fHistPyCM = new TH1F("fHistPyCM","",200,0,2);\r
+ TH1F* fHistPzCM = new TH1F("fHistPzCM","",200,0,2);\r
+ TH1F* fHistKsCM = new TH1F("fHistKsCM","",200,0,2);\r
+ fOutputList->Add(fHistPx);\r
+ fOutputList->Add(fHistPy);\r
+ fOutputList->Add(fHistPz);\r
+ fOutputList->Add(fHistPxCM);\r
+ fOutputList->Add(fHistPyCM);\r
+ fOutputList->Add(fHistPzCM);\r
+ fOutputList->Add(fHistKsCM);\r
+\r
+ TH1F* fHistPOutLCMS = new TH1F("fHistPOutLCMS","",200,0,2);\r
+ TH1F* fHistPSideLCMS = new TH1F("fHistPSideLCMS","",200,0,2);\r
+ TH1F* fHistPLongLCMS = new TH1F("fHistPLongLCMS","",200,0,2);\r
+ fOutputList->Add(fHistPOutLCMS);\r
+ fOutputList->Add(fHistPSideLCMS);\r
+ fOutputList->Add(fHistPLongLCMS);\r
+\r
+ //pair gamma (LCMS to PRF, OSL)\r
+ TH2F* fHistGamma = new TH2F("fHistGamma","Gamma from LCMS to PRF",500,1,5,100,0,1);\r
+ fOutputList->Add(fHistGamma);\r
+\r
//invariant mass distributions\r
TH3F* fHistMass = new TH3F("fHistMass","",kCentBins,.5,kCentBins+.5,50,0.,5.,400,.3,.7);\r
fOutputList->Add(fHistMass);\r
fOutputList->Add(fHistSepNumNeg);\r
TH1F* fHistSepDenNeg = new TH1F("fHistSepDenNeg","",200,0,20);\r
fOutputList->Add(fHistSepDenNeg);\r
- //TH1F* fHistSepNumPosOld = new TH1F("fHistSepNumPosOld","",200,0,20);\r
- //fOutputList->Add(fHistSepNumPosOld);\r
- //TH1F* fHistSepDenPosOld = new TH1F("fHistSepDenPosOld","",200,0,20);\r
- //fOutputList->Add(fHistSepDenPosOld);\r
- //TH1F* fHistSepNumNegOld = new TH1F("fHistSepNumNegOld","",200,0,20);\r
- //fOutputList->Add(fHistSepNumNegOld);\r
- //TH1F* fHistSepDenNegOld = new TH1F("fHistSepDenNegOld","",200,0,20);\r
- //fOutputList->Add(fHistSepDenNegOld);\r
\r
TH2F* fHistSepNumPos2 = new TH2F("fHistSepNumPos2","",100,0,20,100,0,20);\r
TH2F* fHistSepDenPos2 = new TH2F("fHistSepDenPos2","",100,0,20,100,0,20);\r
TH2F* fHistSepNumNeg2 = new TH2F("fHistSepNumNeg2","",100,0,20,100,0,20);\r
TH2F* fHistSepDenNeg2 = new TH2F("fHistSepDenNeg2","",100,0,20,100,0,20);\r
- //TH2F* fHistSepNumPos2Old = new TH2F("fHistSepNumPos2Old","",100,0,20,100,0,20);\r
- //TH2F* fHistSepDenPos2Old = new TH2F("fHistSepDenPos2Old","",100,0,20,100,0,20);\r
- //TH2F* fHistSepNumNeg2Old = new TH2F("fHistSepNumNeg2Old","",100,0,20,100,0,20);\r
- //TH2F* fHistSepDenNeg2Old = new TH2F("fHistSepDenNeg2Old","",100,0,20,100,0,20);\r
fOutputList->Add(fHistSepNumPos2);\r
fOutputList->Add(fHistSepDenPos2);\r
fOutputList->Add(fHistSepNumNeg2);\r
fOutputList->Add(fHistSepDenNeg2);\r
- //fOutputList->Add(fHistSepNumPos2Old);\r
- //fOutputList->Add(fHistSepDenPos2Old);\r
- //fOutputList->Add(fHistSepNumNeg2Old);\r
- //fOutputList->Add(fHistSepDenNeg2Old);\r
+\r
\r
TH2F* fHistSepDPC = new TH2F("fHistSepDPC","",200,-1,1,50,0,10);\r
TH2F* fHistSepDPCBkg = new TH2F("fHistSepDPCBkg","",200,-1,1,50,0,10);\r
fOutputList->Add(fHistSepDPC);\r
- fOutputList->Add(fHistSepDPCBkg);\r
-\r
- TH1F* fHistPx = new TH1F("fHistPx","",200,0,2);\r
- TH1F* fHistPy = new TH1F("fHistPy","",200,0,2);\r
- TH1F* fHistPz = new TH1F("fHistPz","",200,0,2);\r
- TH1F* fHistPxCM = new TH1F("fHistPxCM","",200,0,2);\r
- TH1F* fHistPyCM = new TH1F("fHistPyCM","",200,0,2);\r
- TH1F* fHistPzCM = new TH1F("fHistPzCM","",200,0,2);\r
- TH1F* fHistKsCM = new TH1F("fHistKsCM","",200,0,2);\r
- fOutputList->Add(fHistPx);\r
- fOutputList->Add(fHistPy);\r
- fOutputList->Add(fHistPz);\r
- fOutputList->Add(fHistPxCM);\r
- fOutputList->Add(fHistPyCM);\r
- fOutputList->Add(fHistPzCM);\r
- fOutputList->Add(fHistKsCM);\r
+ fOutputList->Add(fHistSepDPCBkg); \r
\r
/////////Signal Distributions///////////////////\r
\r
\r
\r
//3D out-side-long\r
- TH3F* fHistOSLCentLowKt = new TH3F("fHistOSLCentLowKt","",100,-.5,.5,100,-.5,.5,100,-.5,.5);\r
+ /*TH3F* fHistOSLCentLowKt = new TH3F("fHistOSLCentLowKt","",100,-.5,.5,100,-.5,.5,100,-.5,.5);\r
fOutputList->Add(fHistOSLCentLowKt);\r
TH3F* fHistOSLCentLowKtBkg = new TH3F("fHistOSLCentLowKtBkg","",100,-.5,.5,100,-.5,.5,100,-.5,.5);\r
fOutputList->Add(fHistOSLCentLowKtBkg);\r
TH3F* fHistOSLSemiCentHighKt = new TH3F("fHistOSLSemiCentHighKt","",100,-.5,.5,100,-.5,.5,100,-.5,.5);\r
fOutputList->Add(fHistOSLSemiCentHighKt);\r
TH3F* fHistOSLSemiCentHighKtBkg = new TH3F("fHistOSLSemiCentHighKtBkg","",100,-.5,.5,100,-.5,.5,100,-.5,.5);\r
- fOutputList->Add(fHistOSLSemiCentHighKtBkg);\r
+ fOutputList->Add(fHistOSLSemiCentHighKtBkg);*/\r
+\r
+ //3D out-side-long\r
+ TH3F *fHist3DOSLSignal[10][4];\r
+ TH3F *fHist3DOSLBkg[10][4];\r
+ \r
+ for(int i3D=0;i3D<10;i3D++){\r
+ for(int j3D=0;j3D<4;j3D++){\r
+ TString *histname = new TString("fHist3DOSL");\r
+ *histname += i3D;\r
+ *histname += j3D;\r
+ histname->Append("Signal");\r
+ fHist3DOSLSignal[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);\r
+ fOutputList->Add(fHist3DOSLSignal[i3D][j3D]);\r
+ histname->Replace(12,6,"Bkg");\r
+ fHist3DOSLBkg[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);\r
+ fOutputList->Add(fHist3DOSLBkg[i3D][j3D]);\r
+ }\r
+ }\r
+ \r
+\r
\r
//3D Spherical Harmonics\r
TH3F* fHistSHCentLowKt = new TH3F("fHistSHCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);\r
const float kEtaCut = 0.8; //maximum |pseudorapidity|\r
const float kMinCosAngle = 0.99; //minimum cosine of K0 pointing angle \r
\r
- const float kMinSeparation = 5.0; //minimum daughter (pair) separation\r
+ const float kMinSeparation = fMinSep; //minimum daughter (pair) separation\r
\r
const float kTOFLow = 0.8; //boundary for TOF usage\r
const float kMaxTOFSigmaPion = 3.0; //TOF # of sigmas\r
bool orderswitch = kFALSE;\r
if(tempTrack->Charge() > 0) {pos0or1 = 0; neg0or1 = 1;}\r
else {pos0or1 = 1; neg0or1 = 0; orderswitch = kTRUE;}\r
+ //tempTrack->~AliAODTrack();\r
\r
//load daughter tracks\r
AliAODTrack* prongTrackPos = (AliAODTrack*)v0->GetDaughter(pos0or1);\r
}\r
\r
//K0 cuts\r
- if(v0->Eta() > kEtaCut) continue; \r
- if(v0->CosPointingAngle(primaryVertex) < kMinCosAngle) continue;\r
- if(v0->MassK0Short() < .2 || v0->MassK0Short() > .8) continue;\r
- if(v0->DcaNegToPrimVertex() < kMinDCAPrimaryPion) continue;\r
- if(v0->DcaPosToPrimVertex() < kMinDCAPrimaryPion) continue; \r
- if(v0->DecayLength(primaryVertex) > kMaxDLK0) continue;\r
- if(v0->DecayLength(primaryVertex) < kMinDLK0) continue;\r
- if(v0->DcaV0Daughters() > kMaxDCADaughtersK0) continue;\r
+ if(v0->Eta() > kEtaCut) continue; \r
+ if(v0->CosPointingAngle(primaryVertex) < kMinCosAngle) continue;\r
+ if(v0->MassK0Short() < .2 || v0->MassK0Short() > .8) continue;\r
+ if(v0->DcaNegToPrimVertex() < kMinDCAPrimaryPion) continue;\r
+ if(v0->DcaPosToPrimVertex() < kMinDCAPrimaryPion) continue; \r
+ if(v0->DecayLength(primaryVertex) > kMaxDLK0) continue;\r
+ if(v0->DecayLength(primaryVertex) < kMinDLK0) continue;\r
+ if(v0->DcaV0Daughters() > kMaxDCADaughtersK0) continue;\r
double v0Dca = v0->DcaV0ToPrimVertex();\r
- if(v0Dca > kMaxDCAK0) continue; \r
- if(!goodPiMinus || !goodPiPlus) continue; \r
+ if(v0Dca > kMaxDCAK0) continue; \r
+ if(!goodPiMinus || !goodPiPlus) continue; \r
\r
//EVERYTHING BELOW HERE PASSES SINGLE PARTICLE CUTS, PION PID, and LOOSE MASS CUT\r
\r
//else continue; //removed, Apr 18\r
\r
//Check for shared daughters, using v0 DCA to judge\r
+ bool v0JudgeNew; //true if new v0 beats old\r
tempK0[v0Count].fSkipShared = kFALSE;\r
+ double newV0Pars[3] = {fabs(v0->MassK0Short()-kMassK0Short),v0Dca,v0->DcaV0Daughters()}; //parameters used in merit cut\r
if(fMeritCase){\r
for(int iID = 0; iID<v0Count; iID++){\r
if(tempK0[iID].fSkipShared == kFALSE){ //if old is already skipped, go to next old\r
if(tempK0[iID].fDaughterID1 == prongTrackPos->GetID() || tempK0[iID].fDaughterID2 == prongTrackNeg->GetID()){\r
- if(tempK0[iID].fV0Dca <= v0Dca){ //if old beats new...\r
+ double oldV0Pars[3] = {fabs(tempK0[iID].fMass-kMassK0Short), tempK0[iID].fV0Dca, tempK0[iID].fDDDca}; \r
+ v0JudgeNew = CheckMeritCutWinner(fMeritCutChoice, oldV0Pars, newV0Pars); //true if new wins\r
+ if(!v0JudgeNew){ //if old beats new...\r
if(!tempK0[iID].fK0 && goodK0) continue; //if bad old beats new good, do nothing... \r
else{ //but if bad old beats new bad, or good old beats anything, skip new\r
tempK0[v0Count].fSkipShared = kTRUE; //skip new\r
else{ //if new beats old...\r
if(tempK0[iID].fK0 && !goodK0) continue; //if bad new beats good old, do nothing...\r
else{ //but if bad new beats bad old, or good new beats anything, skip old\r
- tempK0[iID].fSkipShared = kTRUE; //skip old \r
- if(tempK0[iID].fK0) k0Count--; //if good old gets skipped, subtract from number of K0s (new one will be added later, if it succeeds)\r
+ tempK0[iID].fSkipShared = kTRUE; //skip old \r
+ if(tempK0[iID].fK0) k0Count--; //if good old gets skipped, subtract from number of K0s (new one will be added later, if it succeeds)\r
}\r
}\r
}\r
}\r
}\r
if(tempK0[v0Count].fSkipShared) continue; //if new K0 is skipped, don't load; go to next v0\r
- }//if MeritCase \r
- \r
+ }//if MeritCase \r
+ \r
//load parameters into temporary class instance\r
if(v0Count < kMaxNumK0)\r
{\r
//else tempK0[v0Count].fSideLeft = kFALSE;\r
//if(v0->MassK0Short() > .515 && v0->MassK0Short() < .545) tempK0[v0Count].fSideRight = kTRUE;\r
//else tempK0[v0Count].fSideRight = kFALSE;\r
- //if(!goodK0) continue; //no sides, speed up analysis (REDUNDANT RIGHT NOW)\r
+ //if(!goodK0) continue; //no sides, speed up analysis (REDUNDANT RIGHT NOW)\r
\r
tempK0[v0Count].fDaughterID1 = prongTrackPos->GetID();\r
tempK0[v0Count].fDaughterID2 = prongTrackNeg->GetID();\r
\r
//for hists\r
tempK0[v0Count].fDDDca = v0->DcaV0Daughters();\r
- tempK0[v0Count].fDecayLength = v0->DecayLength(primaryVertex);\r
+ tempK0[v0Count].fDecayLength = v0->DecayLength(primaryVertex);\r
tempK0[v0Count].fPosPt = v0->PtProng(pos0or1);\r
tempK0[v0Count].fNegPt = v0->PtProng(neg0or1);\r
tempK0[v0Count].fPosPhi = v0->PhiProng(pos0or1);\r
tempK0[v0Count].fNegPhi = v0->PhiProng(neg0or1);\r
- if(!orderswitch){\r
+ if(!orderswitch){\r
tempK0[v0Count].fPosDca = v0->DcaPosToPrimVertex();\r
tempK0[v0Count].fNegDca = v0->DcaNegToPrimVertex();\r
- }\r
+ }\r
else{\r
tempK0[v0Count].fPosDca = v0->DcaNegToPrimVertex();\r
tempK0[v0Count].fNegDca = v0->DcaPosToPrimVertex();\r
- } \r
-\r
+ } \r
+ \r
//for separation\r
- prongTrackPos->GetCovarianceXYZPxPyPz(tempK0[v0Count].fCovPos);\r
- prongTrackNeg->GetCovarianceXYZPxPyPz(tempK0[v0Count].fCovNeg);\r
- prongTrackPos->GetXYZ(tempK0[v0Count].fXPos);\r
- prongTrackNeg->GetXYZ(tempK0[v0Count].fXNeg);\r
+ GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackPos, bField, tempK0[v0Count].fPosXYZ, vertex);\r
+ GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackNeg, bField, tempK0[v0Count].fNegXYZ, vertex);\r
+ //for DPC\r
prongTrackPos->GetPxPyPz(tempK0[v0Count].fPPos);\r
prongTrackNeg->GetPxPyPz(tempK0[v0Count].fPNeg);\r
- \r
+\r
//if(idCount < 50){\r
// if(goodK0){\r
// idArray[idCount*2] = prongTrackPos->GetID();\r
// Correlations\r
//////////////////////////////////////////////////////////////////////\r
\r
- float px1, py1, pz1, px2, py2, pz2, en1, en2; //single kaon values\r
- float pairPx, pairPy, pairPz; //kaon pair values\r
- float pairP0, pairPt, pairKt, pairMt; //LCMS values for out-side-long\r
- float qinv, q0, qx, qy, qz, qLong, qSide, qOut; //pair q values\r
+ float px1, py1, pz1, px2, py2, pz2; //single kaon values\r
+ float en1, en2; //single kaon values \r
+ //float pt1, pt2; //single kaon values\r
+ float pairPx, pairPy, pairPz, pairP0; //pair momentum values\r
+ float pairPt, pairMt, pairKt; //pair momentum values\r
+ float pairMInv, pairPDotQ;\r
+ float qinv, q0, qx, qy, qz; //pair q values\r
float qLength, thetaSH, thetaSHCos, phiSH; //Spherical Harmonics values\r
- //float pt1, pt2; //single kaon pt\r
float am12, epm, h1, p12, p112, ppx, ppy, ppz, ks; //PRF\r
+ //float qOutLCMS;\r
+ float qOutPRF, qSide, qLong; //relative momentum in LCMS/PRF frame\r
+ float betasq, gamma;\r
+ float p1LCMSOut, p1LCMSSide, p1LCMSLong, en1LCMS;\r
+ float p2LCMSOut, p2LCMSSide, p2LCMSLong, en2LCMS;\r
+\r
\r
for(int i=0; i<(fEvt)->fNumV0s; i++) // Current event V0\r
{\r
if(evnum==0) // Get rid of shared tracks\r
{\r
if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;\r
- if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;\r
- if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;\r
- if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;\r
- }\r
-\r
+ if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;\r
+ if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;\r
+ if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;\r
+ }\r
+ \r
px1 = (fEvt)->fK0Particle[i].fMomentum[0];\r
- px2 = (fEvt+evnum)->fK0Particle[j].fMomentum[0];\r
- py1 = (fEvt)->fK0Particle[i].fMomentum[1];\r
- py2 = (fEvt+evnum)->fK0Particle[j].fMomentum[1];\r
- pz1 = (fEvt)->fK0Particle[i].fMomentum[2];\r
- pz2 = (fEvt+evnum)->fK0Particle[j].fMomentum[2];\r
+ py1 = (fEvt)->fK0Particle[i].fMomentum[1];\r
+ pz1 = (fEvt)->fK0Particle[i].fMomentum[2];\r
//pt1 = (fEvt)->fK0Particle[i].fPt;\r
+ px2 = (fEvt+evnum)->fK0Particle[j].fMomentum[0];\r
+ py2 = (fEvt+evnum)->fK0Particle[j].fMomentum[1];\r
+ pz2 = (fEvt+evnum)->fK0Particle[j].fMomentum[2];\r
//pt2 = (fEvt+evnum)->fK0Particle[j].fPt;\r
+ if(fRandomNumber->Rndm() < .5){ //switch particle order for 3D qout bias\r
+ double tempvar;\r
+ tempvar = px1; px1 = px2; px2 = tempvar;\r
+ tempvar = py1; py1 = py2; py2 = tempvar;\r
+ tempvar = pz1; pz1 = pz2; pz2 = tempvar;\r
+ }\r
\r
- pairPx = px1 + px2;\r
+ en1 = sqrt(pow(px1,2)+pow(py1,2)+pow(pz1,2)+pow(kMassK0Short,2));\r
+ en2 = sqrt(pow(px2,2)+pow(py2,2)+pow(pz2,2)+pow(kMassK0Short,2));\r
\r
- pairPy = py1 + py2;\r
- pairPz = pz1 + pz2;\r
- pairKt = sqrt(pairPx*pairPx + pairPy*pairPy)/2.;\r
+ q0 = en1 - en2;\r
+ qx = px1 - px2;\r
+ qy = py1 - py2;\r
+ qz = pz1 - pz2;\r
+ qinv = sqrt(pow(qx,2) + pow(qy,2) + pow(qz,2) - pow(q0,2));\r
\r
- en1 = sqrt(pow(px1,2)+pow(py1,2)+pow(pz1,2)+pow(kMassK0Short,2));\r
- en2 = sqrt(pow(px2,2)+pow(py2,2)+pow(pz2,2)+pow(kMassK0Short,2));\r
- \r
- qinv = sqrt(pow(px1-px2,2) + pow(py1-py2,2) + pow(pz1-pz2,2) - pow(en1-en2,2));\r
- \r
- //PRF\r
- p12 = sqrt(pow(pairPx,2)+pow(pairPy,2)+pow(pairPz,2)); //pair momentum length\r
- am12 = sqrt(pow(en1+en2,2)-p12*p12); //sqrt(s), |p1+p2| (4vec)\r
- epm = en1+en2+am12; //"energy plus mass"\r
- p112 = px1*pairPx+py1*pairPy+pz1*pairPz; //proj. of p1 on pairP\r
+ pairPx = px1 + px2;\r
+ pairPy = py1 + py2;\r
+ pairPz = pz1 + pz2;\r
+ pairP0 = en1 + en2;\r
+ pairPt = sqrt(pairPx*pairPx + pairPy*pairPy);\r
+ pairKt = pairPt/2.; //used for KT binning\r
+ pairMt = sqrt(pairP0*pairP0 - pairPz*pairPz); //used for LCMS (not plots)\r
+ pairMInv = sqrt(pow(pairP0,2)-pow(pairPx,2)-pow(pairPy,2)-pow(pairPz,2));//used for PRF\r
+ pairPDotQ = pairP0*q0-pairPx*qx-pairPy*qy-pairPz*qz; //used for PRF\r
+\r
+ //PRF (this section will probably be removed in favor of later boosting section)\r
+ p12 = sqrt(pow(pairPx,2)+pow(pairPy,2)+pow(pairPz,2)); //pair momentum length\r
+ am12 = sqrt(pow(en1+en2,2)-p12*p12); //sqrt(s)=|p1+p2|(4vec)\r
+ epm = en1+en2+am12; //"energy plus mass"\r
+ p112 = px1*pairPx+py1*pairPy+pz1*pairPz; //proj. of p1 on pairP\r
+ if(am12 == 0) continue;\r
h1 = (p112/epm - en1)/am12;\r
- ppx = px1+pairPx*h1; //px in PRF\r
- ppy = py1+pairPy*h1; //py in PRF \r
- ppz = pz1+pairPz*h1; //pz in PRF\r
- ks = sqrt(ppx*ppx+ppy*ppy+ppz*ppz); //k*\r
+ ppx = px1+pairPx*h1; //px in PRF\r
+ ppy = py1+pairPy*h1; //py in PRF \r
+ ppz = pz1+pairPz*h1; //pz in PRF\r
+ ks = sqrt(ppx*ppx+ppy*ppy+ppz*ppz); //k*\r
((TH1F*)fOutputList->FindObject("fHistPxCM"))->Fill(ppx);\r
((TH1F*)fOutputList->FindObject("fHistPyCM"))->Fill(ppy);\r
((TH1F*)fOutputList->FindObject("fHistPzCM"))->Fill(ppz);\r
((TH1F*)fOutputList->FindObject("fHistKsCM"))->Fill(ks);\r
\r
-\r
- //out-side-long\r
- pairP0 = en1 + en2;\r
- q0 = en1 - en2;\r
- qx = px1 - px2;\r
- qy = py1 - py2;\r
- qz = pz1 - pz2;\r
- if(fRandomNumber->Rndm() < .5){\r
- //qx = -1*qx; qy = -1*qy; qz = -1*qz;\r
- }\r
- pairPt = pairKt*2.;\r
- pairMt = sqrt(pairP0*pairP0 - pairPz*pairPz);\r
- qLong = (pairP0*qz - pairPz*q0)/pairMt;\r
- qOut = (pairPx*qx + pairPy*qy)/pairPt;\r
- qSide = (pairPx*qy - pairPy*qx)/pairPt;\r
+ //relative momentum in out-side-long for LCMS and PRF\r
+ if(pairMt == 0 || pairPt == 0) continue;\r
+ qLong = (pairP0*qz - pairPz*q0)/pairMt; //same for both frames\r
+ qSide = (pairPx*qy - pairPy*qx)/pairPt; //same for both frames\r
+ //qOutLCMS = (pairPx*qx + pairPy*qy)/pairPt;\r
+ qOutPRF = pairMInv*(pairPx*qx+pairPy*qy)/pairMt/pairPt - pairPt*pairPDotQ/pairMt/pairMInv;\r
+ \r
+ //finding gamma for gamma binning/hists (likely will be removed after tests)\r
+ p1LCMSOut = (pairPx*px1+pairPy*py1)/pairPt;\r
+ p1LCMSSide = (pairPx*py1-pairPy*px1)/pairPt;\r
+ p1LCMSLong = (pairP0*pz1-pairPz*en1)/pairMt;\r
+ p2LCMSOut = (pairPx*px2+pairPy*py2)/pairPt;\r
+ p2LCMSSide = (pairPx*py2-pairPy*px2)/pairPt;\r
+ p2LCMSLong = (pairP0*pz2-pairPz*en2)/pairMt;\r
+ en1LCMS = sqrt(pow(p1LCMSOut,2)+pow(p1LCMSSide,2)+pow(p1LCMSLong,2)+pow(kMassK0Short,2));\r
+ en2LCMS = sqrt(pow(p2LCMSOut,2)+pow(p2LCMSSide,2)+pow(p2LCMSLong,2)+pow(kMassK0Short,2)); \r
+ betasq = pow((p1LCMSOut+p2LCMSOut)/(en1LCMS+en2LCMS),2);\r
+ gamma = 1./sqrt(1-betasq);\r
+ ((TH2F*)fOutputList->FindObject("fHistGamma"))->Fill(gamma,qinv);\r
+ ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p1LCMSOut);\r
+ ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p1LCMSSide);\r
+ ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p1LCMSLong);\r
+ ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p2LCMSOut);\r
+ ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p2LCMSSide);\r
+ ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p2LCMSLong);\r
+ //getting bin numbers and names for 3D histogram\r
+ TString *histname3D = new TString("fHist3DOSL");\r
+ int ktBin;\r
+ if(pairKt < 0.6) ktBin = 0;\r
+ else if(pairKt < 0.8) ktBin = 1;\r
+ else if(pairKt < 1.0) ktBin = 2;\r
+ else ktBin = 3;\r
+ *histname3D += centBin-6; //centBins: [6,15] -> array bins: [0,9]\r
+ *histname3D += ktBin;\r
\r
//Spherical harmonics\r
- qLength = sqrt(qLong*qLong + qSide*qSide + qOut*qOut);\r
+ qLength = sqrt(qLong*qLong + qSide*qSide + qOutPRF*qOutPRF);\r
thetaSHCos = qLong/qLength;\r
thetaSH = acos(thetaSHCos);\r
- phiSH = acos(qOut/(qLength*sin(thetaSH)));\r
-\r
- //SEPARATION STUDIES (two methods are compared here; one will be phased out soon (old way is commented out))\r
- //Both methods take same-sign daughter separation throughout TPC\r
- fPosDaughter1->Set((fEvt)->fK0Particle[i].fXPos, (fEvt)->fK0Particle[i].fPPos, (fEvt)->fK0Particle[i].fCovPos, 1);\r
- fNegDaughter1->Set((fEvt)->fK0Particle[i].fXNeg, (fEvt)->fK0Particle[i].fPNeg, (fEvt)->fK0Particle[i].fCovNeg, -1);\r
- fPosDaughter2->Set((fEvt+evnum)->fK0Particle[j].fXPos, (fEvt+evnum)->fK0Particle[j].fPPos, (fEvt+evnum)->fK0Particle[j].fCovPos, 1);\r
- fNegDaughter2->Set((fEvt+evnum)->fK0Particle[j].fXNeg, (fEvt+evnum)->fK0Particle[j].fPNeg, (fEvt+evnum)->fK0Particle[j].fCovNeg, -1);\r
- \r
- //variables for old method\r
- //double rP1[3]; //positive daughter position (K0 #1)\r
- //double rN1[3]; //negative daughter position (K0 #1)\r
- //double rP2[3]; //positive daughter position (K0 #2)\r
- //double rN2[3]; //negative daughter position (K0 #2)\r
- //float pDiff; //positive daughter separation\r
- //float nDiff; //negative daughter separation\r
- //float pMean = 0; //average separation, positive\r
- //float nMean = 0; //average separation, negative\r
- //float pMin = 9999; //minimum separation, positive\r
- //float nMin = 9999; //minimum separation, negative\r
-\r
- //new method from Hans Beck\r
+ phiSH = acos(qOutPRF/(qLength*sin(thetaSH)));\r
+\r
+ //Finding average separation of daughters throughout TPC - two-track cut\r
float posPositions1[9][3] = {{0}};\r
float negPositions1[9][3] = {{0}};\r
float posPositions2[9][3] = {{0}};\r
float negPositions2[9][3] = {{0}};\r
- GetGlobalPositionAtGlobalRadiiThroughTPC(fPosDaughter1,bField,posPositions1);\r
- GetGlobalPositionAtGlobalRadiiThroughTPC(fPosDaughter2,bField,posPositions2);\r
- GetGlobalPositionAtGlobalRadiiThroughTPC(fNegDaughter1,bField,negPositions1);\r
- GetGlobalPositionAtGlobalRadiiThroughTPC(fNegDaughter2,bField,negPositions2);\r
- float pMean2 = 0;\r
- float nMean2 = 0;\r
-\r
- float pDiff2;\r
- float nDiff2;\r
- float pMin2 = 9999;\r
- float nMin2 = 9999;\r
-\r
- double pCount=0; //counter for number of radial points used (low pT tracks don't go all the way through TPC)\r
- double nCount=0;\r
+ for(int iPos = 0; iPos < 9; iPos++){\r
+ for(int jPos = 0; jPos < 3; jPos++){\r
+ posPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fPosXYZ[iPos][jPos];\r
+ negPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fNegXYZ[iPos][jPos];\r
+ posPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fPosXYZ[iPos][jPos];\r
+ negPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fNegXYZ[iPos][jPos];\r
+ }\r
+ } \r
+ float pMean = 0.; //average separation for positive daughters\r
+ float nMean = 0.; //average separation for negative daughters\r
+ float pDiff; \r
+ float nDiff;\r
+ float pMin = 9999.; //minimum separation (updates) - pos\r
+ float nMin = 9999.; //minimum separation (updates) - neg\r
+ double pCount=0; //counter for number of points used - pos\r
+ double nCount=0; //counter for number of points used - neg\r
for(int ss=0;ss<9;ss++){\r
- \r
if(posPositions1[ss][0] != -9999 && posPositions2[ss][0] != -9999){ \r
pCount++;\r
- //fPosDaughter1->GetXYZAt(85+(ss*20), bField, rP1);\r
- //fPosDaughter2->GetXYZAt(85+(ss*20), bField, rP2);\r
- //pDiff = sqrt(pow(rP1[0]-rP2[0],2)+pow(rP1[1]-rP2[1],2)+pow(rP1[2]-rP2[2],2));\r
- pDiff2 = sqrt(pow(posPositions1[ss][0]-posPositions2[ss][0],2)+pow(posPositions1[ss][1]-posPositions2[ss][1],2)+pow(posPositions1[ss][2]-posPositions2[ss][2],2));\r
- //pMean = pMean + pDiff;\r
- pMean2 = pMean2 + pDiff2;\r
- //if(pDiff < pMin) pMin = pDiff;\r
- if(pDiff2 < pMin2) pMin2 = pDiff2;\r
+ pDiff = sqrt(pow(posPositions1[ss][0]-posPositions2[ss][0],2)+pow(posPositions1[ss][1]-posPositions2[ss][1],2)+pow(posPositions1[ss][2]-posPositions2[ss][2],2));\r
+ pMean = pMean + pDiff;\r
+ if(pDiff < pMin) pMin = pDiff;\r
}\r
-\r
if(negPositions1[ss][0] != -9999 && negPositions1[ss][0] != -9999){\r
nCount++;\r
- //fNegDaughter1->GetXYZAt(85+(ss*20), bField, rN1);\r
- //fNegDaughter2->GetXYZAt(85+(ss*20), bField, rN2);\r
- //nDiff = sqrt(pow(rN1[0]-rN2[0],2)+pow(rN1[1]-rN2[1],2)+pow(rN1[2]-rN2[2],2));\r
- nDiff2 = sqrt(pow(negPositions1[ss][0]-negPositions2[ss][0],2)+pow(negPositions1[ss][1]-negPositions2[ss][1],2)+pow(negPositions1[ss][2]-negPositions2[ss][2],2)); \r
- //nMean = nMean + nDiff;\r
- nMean2 = nMean2 + nDiff2;\r
- //if(nDiff < nMin) nMin = nDiff;\r
- if(nDiff2 < nMin2) nMin2 = nDiff2;\r
+ nDiff = sqrt(pow(negPositions1[ss][0]-negPositions2[ss][0],2)+pow(negPositions1[ss][1]-negPositions2[ss][1],2)+pow(negPositions1[ss][2]-negPositions2[ss][2],2)); \r
+ nMean = nMean + nDiff;\r
+ if(nDiff < nMin) nMin = nDiff;\r
}\r
}\r
- //pMean = pMean/pCount;\r
- //nMean = nMean/nCount;\r
- pMean2 = pMean2/pCount;\r
- nMean2 = nMean2/nCount; \r
+ pMean = pMean/pCount; \r
+ nMean = nMean/nCount; \r
\r
if(evnum==0){ \r
- ((TH1F*)fOutputList->FindObject("fHistSepNumPos"))->Fill(pMean2); \r
- ((TH1F*)fOutputList->FindObject("fHistSepNumNeg"))->Fill(nMean2);\r
- //((TH1F*)fOutputList->FindObject("fHistSepNumPosOld"))->Fill(pMean);\r
- //((TH1F*)fOutputList->FindObject("fHistSepNumNegOld"))->Fill(nMean);\r
- ((TH2F*)fOutputList->FindObject("fHistSepNumPos2"))->Fill(pMean2,pMin2);\r
- ((TH2F*)fOutputList->FindObject("fHistSepNumNeg2"))->Fill(nMean2,nMin2);\r
- //((TH2F*)fOutputList->FindObject("fHistSepNumPos2Old"))->Fill(pMean,pMin);\r
- //((TH2F*)fOutputList->FindObject("fHistSepNumNeg2Old"))->Fill(nMean,nMin);\r
+ ((TH1F*)fOutputList->FindObject("fHistSepNumPos"))->Fill(pMean); \r
+ ((TH1F*)fOutputList->FindObject("fHistSepNumNeg"))->Fill(nMean);\r
+ ((TH2F*)fOutputList->FindObject("fHistSepNumPos2"))->Fill(pMean,pMin);\r
+ ((TH2F*)fOutputList->FindObject("fHistSepNumNeg2"))->Fill(nMean,nMin);\r
}\r
else{\r
- ((TH1F*)fOutputList->FindObject("fHistSepDenPos"))->Fill(pMean2); \r
- ((TH1F*)fOutputList->FindObject("fHistSepDenNeg"))->Fill(nMean2);\r
- //((TH1F*)fOutputList->FindObject("fHistSepDenPosOld"))->Fill(pMean);\r
- //((TH1F*)fOutputList->FindObject("fHistSepDenNegOld"))->Fill(nMean);\r
- ((TH2F*)fOutputList->FindObject("fHistSepDenPos2"))->Fill(pMean2,pMin2);\r
- ((TH2F*)fOutputList->FindObject("fHistSepDenNeg2"))->Fill(nMean2,nMin2);\r
- //((TH2F*)fOutputList->FindObject("fHistSepDenPos2Old"))->Fill(pMean,pMin);\r
- //((TH2F*)fOutputList->FindObject("fHistSepDenNeg2Old"))->Fill(nMean,nMin);\r
- }\r
+ ((TH1F*)fOutputList->FindObject("fHistSepDenPos"))->Fill(pMean); \r
+ ((TH1F*)fOutputList->FindObject("fHistSepDenNeg"))->Fill(nMean);\r
+ ((TH2F*)fOutputList->FindObject("fHistSepDenPos2"))->Fill(pMean,pMin);\r
+ ((TH2F*)fOutputList->FindObject("fHistSepDenNeg2"))->Fill(nMean,nMin);\r
+ }\r
\r
//Decay plane coincidence\r
//daughter momenta\r
float crosslength2 = sqrt(pow(cross2[0],2)+pow(cross2[1],2)+pow(cross2[2],2));\r
float dpc = (cross1[0]*cross2[0]+cross1[1]*cross2[1]+cross1[2]*cross2[2])/(crosslength1*crosslength2);\r
\r
- if(evnum==0)((TH2F*)fOutputList->FindObject("fHistSepDPC"))->Fill(dpc,pMean2);\r
- else ((TH2F*)fOutputList->FindObject("fHistSepDPCBkg"))->Fill(dpc,pMean2);\r
+ if(evnum==0)((TH2F*)fOutputList->FindObject("fHistSepDPC"))->Fill(dpc,pMean);\r
+ else ((TH2F*)fOutputList->FindObject("fHistSepDPCBkg"))->Fill(dpc,pMean);\r
\r
- if(pMean2 < kMinSeparation || nMean2 < kMinSeparation) continue; //using the "new" method (ala Hans)\r
+ if(pMean < kMinSeparation || nMean < kMinSeparation) continue; //using the "new" method (ala Hans)\r
//end separation studies\r
\r
//Fill Histograms\r
bool center1K0 = kFALSE; //accepted mass K0\r
- bool center2K0 = kFALSE;\r
+ bool center2K0 = kFALSE;\r
if((fEvt)->fK0Particle[i].fK0) center1K0=kTRUE;\r
if((fEvt+evnum)->fK0Particle[j].fK0) center2K0=kTRUE;\r
//bool CL1 = kFALSE;\r
//else if((fEvt)->fK0Particle[i].fSideRight) SideRight1 = kTRUE;\r
//if((fEvt+evnum)->fK0Particle[j].fSideLeft) SideLeft2 = kTRUE;\r
//else if((fEvt+evnum)->fK0Particle[j].fSideRight) SideRight2 = kTRUE;\r
-\r
- //for daughter sharing studies - REMOVED - NOW I CUT SHARED DAUGHTERS AT THE V0 LEVEL\r
- //bool splitK0sides = kFALSE;\r
- //bool splitK0centers = kFALSE;\r
- //int posD1ID = (fEvt)->fK0Particle[i].fDaughterID1;\r
- //int negD1ID = (fEvt)->fK0Particle[i].fDaughterID2;\r
- //int posD2ID = (fEvt+evnum)->fK0Particle[j].fDaughterID1;\r
- //int negD2ID = (fEvt+evnum)->fK0Particle[j].fDaughterID2;\r
- //if(evnum == 0){\r
- // //centers\r
- // if(center1K0 && center2K0){\r
- // for(int iID=0;iID<idCount;iID++){\r
- // if(posD1ID == idArray[iID*2] && negD2ID == idArray[iID*2+1]){\r
- // ((TH3F*)fOutputList->FindObject("fHistSplitK0Centers"))->Fill(centBin+1, pairKt, qinv);\r
- // splitK0centers = kTRUE;}\r
- // else if(negD1ID == idArray[iID*2+1] && posD2ID == idArray[iID*2]){\r
- // ((TH3F*)fOutputList->FindObject("fHistSplitK0Centers"))->Fill(centBin+1, pairKt, qinv);\r
- // splitK0centers = kTRUE;}\r
- // }}\r
- // //sides\r
- // else if((SideLeft1 || SideRight1) && (SideLeft2 || SideRight2)){\r
- // for(int iID=0;iID<idCount;iID++){\r
- // if(posD1ID == idArray[iID*2] && negD2ID == idArray[iID*2+1]){\r
- // ((TH3F*)fOutputList->FindObject("fHistSplitK0Sides"))->Fill(centBin+1, pairKt, qinv);\r
- // splitK0sides = kTRUE;}\r
- // else if(negD1ID == idArray[iID*2+1] && posD2ID == idArray[iID*2]){\r
- // ((TH3F*)fOutputList->FindObject("fHistSplitK0Sides"))->Fill(centBin+1, pairKt, qinv);\r
- // splitK0sides = kTRUE;}\r
- // }}\r
- //}//end of daughter sharing section\r
\r
\r
if(evnum==0) //Same Event\r
//else ((TH3F*)fOutputList->FindObject("fHistCRCRSignal"))->Fill(centBin+1, pairKt, qinv);\r
\r
//3D\r
- if(pairKt < 1.0){\r
+ if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){\r
+ histname3D->Append("Signal");\r
+ ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);\r
+ }\r
+ \r
+ /*if(pairKt < 1.0){\r
if(centBin > 13){\r
- ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKt"))->Fill(qOut,qSide,qLong);\r
+ ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKt"))->Fill(qOutPRF,qSide,qLong);\r
((TH3F*)fOutputList->FindObject("fHistSHCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}\r
else if(centBin > 9){\r
- ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKt"))->Fill(qOut,qSide,qLong);\r
+ ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKt"))->Fill(qOutPRF,qSide,qLong);\r
((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}} \r
else if(pairKt < 2.0){\r
if(centBin > 13){\r
- ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKt"))->Fill(qOut,qSide,qLong);\r
+ ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKt"))->Fill(qOutPRF,qSide,qLong);\r
((TH3F*)fOutputList->FindObject("fHistSHCentHighKt"))->Fill(qLength,thetaSHCos, phiSH);}\r
else if(centBin > 9){\r
- ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKt"))->Fill(qOut,qSide,qLong);\r
- ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKt"))->Fill(qLength, thetaSHCos, phiSH);}}\r
- }\r
+ ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKt"))->Fill(qOutPRF,qSide,qLong);\r
+\r
+ ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKt"))->Fill(qLength, thetaSHCos, phiSH);}}*/ \r
+\r
+ }//centercenter\r
\r
//side-side correlations\r
//if(!splitK0sides){\r
//else ((TH3F*)fOutputList->FindObject("fHistCRCRBkg"))->Fill(centBin+1, pairKt, qinv);\r
\r
//3D\r
- if(pairKt < 1.0){\r
+ if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){\r
+ histname3D->Replace(12,6,"Bkg");\r
+ ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);\r
+ }\r
+ /*if(pairKt < 1.0){\r
if(centBin > 13){\r
- ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKtBkg"))->Fill(qOut,qSide,qLong);\r
+ ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);\r
((TH3F*)fOutputList->FindObject("fHistSHCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}\r
else if(centBin > 9){\r
- ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKtBkg"))->Fill(qOut,qSide,qLong);\r
+ ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);\r
((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}}\r
else if(pairKt < 2.0){\r
if(centBin > 13){\r
- ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKtBkg"))->Fill(qOut,qSide,qLong);\r
+ ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);\r
((TH3F*)fOutputList->FindObject("fHistSHCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}\r
else if(centBin > 9){\r
- ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKtBkg"))->Fill(qOut,qSide,qLong);\r
- ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}}\r
+ ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);\r
+ ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}}*/\r
}\r
\r
//side-side correlations\r
}\r
\r
//_________________________________________________________________________\r
-void AliFemtoK0Analysis::GetGlobalPositionAtGlobalRadiiThroughTPC(const AliESDtrack *track, const Float_t bfield, Float_t globalPositionsAtRadii[9][3]){\r
+void AliFemtoK0Analysis::GetGlobalPositionAtGlobalRadiiThroughTPC(const AliAODTrack *track, const Float_t bfield, Float_t globalPositionsAtRadii[9][3], double PrimaryVertex[3]){\r
// Gets the global position of the track at nine different radii in the TPC\r
// track is the track you want to propagate\r
// bfield is the magnetic field of your event\r
globalPositionsAtRadii[iR][0]=xyz[0];\r
globalPositionsAtRadii[iR][1]=xyz[1];\r
globalPositionsAtRadii[iR][2]=xyz[2];\r
+ //subtract primary vertex, "zero" track for correct mixed-event comparison\r
+ globalPositionsAtRadii[iR][0] -= PrimaryVertex[0];\r
+ globalPositionsAtRadii[iR][1] -= PrimaryVertex[1];\r
+ globalPositionsAtRadii[iR][2] -= PrimaryVertex[2];\r
+\r
// Indicate we want the next radius \r
iR+=1;\r
}\r
}\r
}\r
\r
+bool AliFemtoK0Analysis::CheckMeritCutWinner(int cutChoice, double oldPars[3], double newPars[3]){\r
+ //performs "merit cut" judgement check on v0s with shared daughters, using one of three criteria.\r
+ //if cutChoice = 4, it uses all three criteria, needed 2 of 3 'points'\r
+\r
+ bool newV0Wins = kFALSE;\r
+ double pardiff[3] = {newPars[0]-oldPars[0],\r
+ newPars[1]-oldPars[1],\r
+ newPars[2]-oldPars[2]};\r
+ if(cutChoice > 0 && cutChoice < 4){\r
+ if(pardiff[cutChoice] <= 0.) newV0Wins = kTRUE;\r
+ }\r
+ else if(cutChoice == 4){\r
+ int newWinCount = 0;\r
+ for(int i=0;i<3;i++){if(pardiff[i+1] <= 0) newWinCount++;}\r
+ if(newWinCount > 1) newV0Wins = kTRUE; \r
+ }\r
+ else{};\r
+ return newV0Wins;\r
+}\r
+\r
+\r
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff