K0s code update (Matt Steinpreis)

diff --git a/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx b/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx
index f587e34605cb2fa4ae1bbc51ddd08b3c98ef6924..a00abdb7dd1ae16877914172c6bf162a7f5e50b0 100644 (file)
--- a/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx
+++ b/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx
@@ -40,6 +40,7 @@
 //     - 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
 ////////////////////////////////////////////////////////////////////////////////\r
 \r
 \r
@@ -255,15 +256,16 @@ void AliFemtoK0Analysis::UserCreateOutputObjects()
   fOutputList->Add(fHistPzCM);\r
   fOutputList->Add(fHistKsCM);\r
 \r
-  TH1F* fHistPtLCMS    = new TH1F("fHistPtLCMS","",300,0,3);\r
-  fOutputList->Add(fHistPtLCMS);\r
-  TH3F *fHistPLCMS = new TH3F("fHistPLCMS","",200,0,2,200,0,2,200,0,2);\r
-  fOutputList->Add(fHistPLCMS);\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
-  TH3F* fHistKtGammaQinv = new TH3F("fHistKtGammaQinv","",200,0,2,500,1,5,100,0,1);\r
-  fOutputList->Add(fHistKtGammaQinv);\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
@@ -396,7 +398,7 @@ void AliFemtoK0Analysis::Exec(Option_t *)
 {\r
   // Main loop\r
   // Called for each event\r
-  cout<<"===========  Event # "<<fEventCount+1<<"  ==========="<<endl;\r
+  //cout<<"===========  Event # "<<fEventCount+1<<"  ==========="<<endl;\r
   fEventCount++;\r
   fAOD = dynamic_cast<AliAODEvent*> (InputEvent());\r
   if (!fAOD) {Printf("ERROR: fAOD not available"); return;}\r
@@ -611,17 +613,17 @@ void AliFemtoK0Analysis::Exec(Option_t *)
     }\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->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->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
@@ -674,7 +676,7 @@ void AliFemtoK0Analysis::Exec(Option_t *)
     //load parameters into temporary class instance\r
     if(v0Count < kMaxNumK0)\r
     {\r
-               if(goodK0){\r
+       if(goodK0){\r
          tempK0[v0Count].fK0 = kTRUE;\r
          k0Count++;\r
         }\r
@@ -684,7 +686,7 @@ void AliFemtoK0Analysis::Exec(Option_t *)
         //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
@@ -697,23 +699,23 @@ void AliFemtoK0Analysis::Exec(Option_t *)
 \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
         //for separation\r
         GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackPos, bField, tempK0[v0Count].fPosXYZ, vertex);\r
-               GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackNeg, bField, tempK0[v0Count].fNegXYZ, 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
@@ -764,7 +766,8 @@ void AliFemtoK0Analysis::Exec(Option_t *)
   //////////////////////////////////////////////////////////////////////\r
 \r
   float px1, py1, pz1, px2, py2, pz2;                  //single kaon values\r
-  float en1, en2;//, pt1, pt2;                                 //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
@@ -774,8 +777,8 @@ void AliFemtoK0Analysis::Exec(Option_t *)
   //float qOutLCMS;\r
   float qOutPRF, qSide, qLong;                         //relative momentum in LCMS/PRF frame\r
   float betasq, gamma;\r
-  float p1LCMSOut, p1LCMSSide, p1LCMSLong, p1LCMSPt, en1LCMS;\r
-  float p2LCMSOut, p2LCMSSide, p2LCMSLong, p2LCMSPt, en2LCMS;\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
@@ -783,18 +786,18 @@ void AliFemtoK0Analysis::Exec(Option_t *)
     //single particle histograms (done here to avoid "skipped" v0s\r
      ((TH1F*)fOutputList->FindObject("fHistDCADaughters"))     ->Fill((fEvt)->fK0Particle[i].fDDDca);\r
      ((TH1F*)fOutputList->FindObject("fHistDecayLengthK0"))    ->Fill((fEvt)->fK0Particle[i].fDecayLength);\r
-     ((TH1F*)fOutputList->FindObject("fHistDCAK0"))                    ->Fill((fEvt)->fK0Particle[i].fV0Dca);\r
+     ((TH1F*)fOutputList->FindObject("fHistDCAK0"))            ->Fill((fEvt)->fK0Particle[i].fV0Dca);\r
      ((TH1F*)fOutputList->FindObject("fHistDCAPiMinus"))       ->Fill((fEvt)->fK0Particle[i].fNegDca);\r
      ((TH1F*)fOutputList->FindObject("fHistDCAPiPlus"))                ->Fill((fEvt)->fK0Particle[i].fPosDca);\r
-     ((TH2F*)fOutputList->FindObject("fHistPtK0"))                     ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPt);\r
+     ((TH2F*)fOutputList->FindObject("fHistPtK0"))             ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPt);\r
      ((TH2F*)fOutputList->FindObject("fHistK0PiPlusPt"))       ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPosPt);\r
      ((TH2F*)fOutputList->FindObject("fHistK0PiMinusPt"))      ->Fill(centBin+1, (fEvt)->fK0Particle[i].fNegPt);\r
      ((TH1F*)fOutputList->FindObject("fHistDaughterPhi"))      ->Fill((fEvt)->fK0Particle[i].fPosPhi);\r
      ((TH1F*)fOutputList->FindObject("fHistDaughterPhi"))      ->Fill((fEvt)->fK0Particle[i].fNegPhi);\r
      \r
-     ((TH1F*)fOutputList->FindObject("fHistPx"))->Fill((fEvt)->fK0Particle[i].fMomentum[0]);\r
-     ((TH1F*)fOutputList->FindObject("fHistPy"))->Fill((fEvt)->fK0Particle[i].fMomentum[1]);\r
-     ((TH1F*)fOutputList->FindObject("fHistPz"))->Fill((fEvt)->fK0Particle[i].fMomentum[2]);\r
+     ((TH1F*)fOutputList->FindObject("fHistPx"))               ->Fill((fEvt)->fK0Particle[i].fMomentum[0]);\r
+     ((TH1F*)fOutputList->FindObject("fHistPy"))               ->Fill((fEvt)->fK0Particle[i].fMomentum[1]);\r
+     ((TH1F*)fOutputList->FindObject("fHistPz"))               ->Fill((fEvt)->fK0Particle[i].fMomentum[2]);\r
 \r
     for(int evnum=0; evnum<kEventsToMix+1; evnum++)// Event buffer loop: evnum=0 is the current event, all other evnum's are past events\r
     {\r
@@ -810,15 +813,22 @@ void AliFemtoK0Analysis::Exec(Option_t *)
              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
+            \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
+        //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
-               //pt1 = (fEvt)->fK0Particle[i].fPt;\r
-               //pt2 = (fEvt+evnum)->fK0Particle[j].fPt;\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
                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
@@ -855,9 +865,6 @@ void AliFemtoK0Analysis::Exec(Option_t *)
         ((TH1F*)fOutputList->FindObject("fHistKsCM"))->Fill(ks);\r
         \r
         //relative momentum in out-side-long for LCMS and PRF\r
-        //if(fRandomNumber->Rndm() < .5){\r
-                //qx = -1*qx; qy = -1*qy; qz = -1*qz;  //randomizing signs, not sure if needed\r
-               //}\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
@@ -865,23 +872,23 @@ void AliFemtoK0Analysis::Exec(Option_t *)
                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
-               p1LCMSPt        = sqrt(pow(p1LCMSOut,2)+pow(p1LCMSSide,2));\r
-               p2LCMSOut       = (pairPx*px2+pairPy*py2)/pairPt;\r
-               p2LCMSSide      = (pairPx*py2-pairPy*px2)/pairPt;\r
-               p2LCMSLong      = (pairP0*pz2-pairPz*en2)/pairMt;\r
-               p2LCMSPt        = sqrt(pow(p2LCMSOut,2)+pow(p2LCMSSide,2));\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
-               ((TH3F*)fOutputList->FindObject("fHistKtGammaQinv"))->Fill(pairKt,gamma,qinv);\r
-               ((TH1F*)fOutputList->FindObject("fHistPtLCMS"))->Fill(p1LCMSPt);\r
-               ((TH1F*)fOutputList->FindObject("fHistPtLCMS"))->Fill(p2LCMSPt);\r
-               ((TH3F*)fOutputList->FindObject("fHistPLCMS"))->Fill(p1LCMSOut,p1LCMSSide,p1LCMSLong);\r
-               ((TH3F*)fOutputList->FindObject("fHistPLCMS"))->Fill(p2LCMSOut,p2LCMSSide,p2LCMSLong);\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
 \r
         //Spherical harmonics\r
         qLength = sqrt(qLong*qLong + qSide*qSide + qOutPRF*qOutPRF);\r
@@ -975,7 +982,7 @@ void AliFemtoK0Analysis::Exec(Option_t *)
 \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
@@ -997,36 +1004,6 @@ void AliFemtoK0Analysis::Exec(Option_t *)
         //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