K0s Analysis code update (Matt Steinpreis)

diff --git a/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx b/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx
index ee3041d130433ec9fcd7c9ce7c1569ae2dd2a663..b42516f0918157264f6dc47c74dfc2d96713fd5c 100644 (file)
--- a/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx
+++ b/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx
@@ -19,10 +19,20 @@
 //  are reconstructed using the AliAODv0 class.  \r
 //\r
 //  authors: Matthew Steinpreis (matthew.steinpreis@cern.ch)\r
+//   \r
+//  Major changes:\r
+//     - TOF mismatch function calls changed (3/28/13)\r
+//     - added minimum decay length cut (3/28/13)\r
 //\r
+//  Minor changes:\r
+//     - K0 multiplicity histogram now filled with "unskippedCount" instead\r
+//       of k0Count (which included skipped k0s with shared daughters) (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
 ////////////////////////////////////////////////////////////////////////////////\r
 \r
 \r
+\r
 #include <iostream>\r
 #include <math.h>\r
 #include "TMath.h"\r
@@ -66,6 +76,7 @@ AliAnalysisTaskSE(),
   fFieldPos(kTRUE),\r
   fOnlineCase(kTRUE),\r
   fMeritCase(kTRUE),\r
+  fMinDecayLength(0.0),\r
   fEventCount(0),\r
   fEC(0x0),\r
   fEvt(0X0),\r
@@ -74,7 +85,6 @@ AliAnalysisTaskSE(),
   fAOD(0x0),\r
   fOutputList(0x0),\r
   fPidAOD(0x0),\r
-  fPidESD(0x0),\r
   fPosDaughter1(0x0),  \r
   fPosDaughter2(0x0),\r
   fNegDaughter1(0x0),\r
@@ -82,11 +92,12 @@ AliAnalysisTaskSE(),
 {\r
 }\r
 //________________________________________________________________________\r
-AliFemtoK0Analysis::AliFemtoK0Analysis(const char *name, bool FieldPositive, bool OnlineCase, bool MeritCase) \r
+AliFemtoK0Analysis::AliFemtoK0Analysis(const char *name, bool FieldPositive, bool OnlineCase, bool MeritCase, float MinDL) \r
 : AliAnalysisTaskSE(name), \r
   fFieldPos(FieldPositive),\r
   fOnlineCase(OnlineCase),\r
   fMeritCase(MeritCase),\r
+  fMinDecayLength(MinDL),\r
   fEventCount(0),\r
   fEC(0x0),\r
   fEvt(0X0),\r
@@ -95,7 +106,6 @@ AliFemtoK0Analysis::AliFemtoK0Analysis(const char *name, bool FieldPositive, boo
   fAOD(0x0),\r
   fOutputList(0x0),\r
   fPidAOD(0x0),\r
-  fPidESD(0x0),\r
   fPosDaughter1(0x0),  \r
   fPosDaughter2(0x0),\r
   fNegDaughter1(0x0),\r
@@ -105,6 +115,8 @@ AliFemtoK0Analysis::AliFemtoK0Analysis(const char *name, bool FieldPositive, boo
   fFieldPos    = FieldPositive;\r
   fOnlineCase  = OnlineCase;\r
   fMeritCase   = MeritCase;\r
+  fMinDecayLength = MinDL;\r
+\r
   // Define output slots here \r
   // Output slot #1\r
   DefineOutput(1, TList::Class());\r
@@ -116,6 +128,7 @@ AliFemtoK0Analysis::AliFemtoK0Analysis(const AliFemtoK0Analysis &obj)
   fFieldPos(obj.fFieldPos),\r
   fOnlineCase(obj.fOnlineCase),\r
   fMeritCase(obj.fMeritCase),\r
+  fMinDecayLength(obj.fMinDecayLength),\r
   fEventCount(obj.fEventCount),\r
   fEC(obj.fEC),\r
   fEvt(obj.fEvt),\r
@@ -124,7 +137,6 @@ AliFemtoK0Analysis::AliFemtoK0Analysis(const AliFemtoK0Analysis &obj)
   fAOD(obj.fAOD),\r
   fOutputList(obj.fOutputList),\r
   fPidAOD(obj.fPidAOD),\r
-  fPidESD(obj.fPidESD),\r
   fPosDaughter1(obj.fPosDaughter1),  \r
   fPosDaughter2(obj.fPosDaughter2),\r
   fNegDaughter1(obj.fNegDaughter1),\r
@@ -140,6 +152,7 @@ AliFemtoK0Analysis &AliFemtoK0Analysis::operator=(const AliFemtoK0Analysis &obj)
  fFieldPos     = obj.fFieldPos;\r
  fOnlineCase   = obj.fOnlineCase;\r
  fMeritCase    = obj.fMeritCase;\r
+ fMinDecayLength= obj.fMinDecayLength;\r
  fEventCount   = obj.fEventCount;\r
  fEC           = obj.fEC;\r
  fEvt          = obj.fEvt;\r
@@ -148,7 +161,6 @@ AliFemtoK0Analysis &AliFemtoK0Analysis::operator=(const AliFemtoK0Analysis &obj)
  fAOD          = obj.fAOD;\r
  fOutputList   = obj.fOutputList;\r
  fPidAOD       = obj.fPidAOD;\r
- fPidESD       = obj.fPidESD;\r
  fPosDaughter1         = obj.fPosDaughter1;  \r
  fPosDaughter2         = obj.fPosDaughter2;\r
  fNegDaughter1         = obj.fNegDaughter1;\r
@@ -167,7 +179,6 @@ AliFemtoK0Analysis::~AliFemtoK0Analysis()
   if(fAOD) delete fAOD;\r
   if(fOutputList) delete fOutputList;\r
   if(fPidAOD) delete fPidAOD;\r
-  if(fPidESD) delete fPidESD;\r
   if(fPosDaughter1) delete fPosDaughter1;\r
   if(fPosDaughter2) delete fPosDaughter2;\r
   if(fNegDaughter1) delete fNegDaughter1;\r
@@ -191,11 +202,8 @@ void AliFemtoK0Analysis::MyInit()
     }\r
   }\r
 \r
-  //fPidAOD = new AliAODpidUtil();\r
   AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
-  fPidAOD = aodH->GetAODpidUtil();\r
-  //fPidAOD = aodH->GetPIDResponse();\r
-  fPidESD = new AliESDpid();\r
+  fPidAOD = aodH->GetPIDResponse();\r
 \r
   fPosDaughter1 = new AliESDtrack();\r
   fPosDaughter2 = new AliESDtrack();\r
@@ -300,6 +308,21 @@ void AliFemtoK0Analysis::UserCreateOutputObjects()
   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
+\r
 /////////Signal Distributions///////////////////\r
 \r
   //1D Q invariant\r
@@ -492,6 +515,7 @@ void AliFemtoK0Analysis::Exec(Option_t *)
   const float kMaxDCADaughtersK0 = 0.3;            //maximum dca of pions to each other - 3D\r
   const float kMaxDCAK0 = 0.3;                     //maximum dca of K0 to primary\r
   const float kMaxDLK0 = 30.0;                     //maximum decay length of K0\r
+  const float kMinDLK0 = fMinDecayLength;         //minimum decay length of K0\r
   const float kEtaCut = 0.8;                       //maximum |pseudorapidity|\r
   const float kMinCosAngle = 0.99;                 //minimum cosine of K0 pointing angle     \r
   \r
@@ -571,24 +595,37 @@ void AliFemtoK0Analysis::Exec(Option_t *)
     if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackNeg,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiMinus = kTRUE;\r
    \r
     //Positive daughter identification TOF\r
+    //float probMis;\r
+    //AliPIDResponse::EDetPidStatus statusPosTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackPos);\r
     double Ppos = v0->PProng(pos0or1);\r
     if(Ppos > kTOFLow) //PiPlus\r
     {\r
      if( (statusPos&AliESDtrack::kTOFpid)!=0 && (statusPos&AliESDtrack::kTIME)!=0 && (statusPos&AliESDtrack::kTOFout)!=0 && (statusPos&AliESDtrack::kTOFmismatch)<=0)\r
+     //if(AliPIDResponse::kDetPidOk == statusPosTOF)\r
      {\r
-      if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackPos,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiPlus = kTRUE;\r
-      else goodPiPlus = kFALSE;\r
-     }  \r
+      //probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);\r
+      //if(probMis < 0.01) //avoid TOF-TPC mismatch\r
+      //{\r
+       if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackPos,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiPlus = kTRUE;\r
+       else goodPiPlus = kFALSE;\r
+      //}  \r
+     }\r
     }\r
     //Negative daughter identification TOF\r
+    //AliPIDResponse::EDetPidStatus statusNegTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackNeg);\r
     double Pneg = v0->PProng(neg0or1);\r
     if(Pneg > kTOFLow) //PiMinus\r
     {\r
      if( (statusNeg&AliESDtrack::kTOFpid)!=0 && (statusNeg&AliESDtrack::kTIME)!=0 && (statusNeg&AliESDtrack::kTOFout)!=0 && (statusNeg&AliESDtrack::kTOFmismatch)<=0)\r
+     //if(AliPIDResponse::kDetPidOk == statusNegTOF)\r
      {\r
-      if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackNeg,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiMinus = kTRUE;\r
-      else goodPiMinus = kFALSE;\r
-      }\r
+      //probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);\r
+      //if(probMis < 0.01) //avoid TOF-TPC mismatch\r
+      //{\r
+       if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackNeg,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiMinus = kTRUE;\r
+       else goodPiMinus = kFALSE;\r
+      //}\r
+     }\r
     }\r
     \r
     //K0 cuts\r
@@ -598,6 +635,7 @@ void AliFemtoK0Analysis::Exec(Option_t *)
     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
@@ -649,13 +687,13 @@ void AliFemtoK0Analysis::Exec(Option_t *)
     {\r
        tempK0[v0Count].fK0 = kTRUE;\r
         //else tempK0[v0Count].fK0 = kFALSE;  //in case I include v0s that arent "good" K0s \r
-        k0Count++;             \r
+        k0Count++; //same as v0count right now (continue if not in mass range, above)\r
 \r
         //if(v0->MassK0Short() > .45 && v0->MassK0Short() < .48) tempK0[v0Count].fSideLeft = kTRUE;\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 \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
@@ -721,7 +759,7 @@ void AliFemtoK0Analysis::Exec(Option_t *)
   //Printf("Number of K0s: %d", k0Count);\r
   tempK0->~AliFemtoK0Particle();\r
 \r
-  ((TH1F*)fOutputList->FindObject("fHistMultK0"))->Fill(k0Count);\r
+  ((TH1F*)fOutputList->FindObject("fHistMultK0"))->Fill(unskippedCount);       // changed 3/25, used to be "k0Count"\r
 \r
   //Printf("Reconstruction Finished. Starting pair studies.");\r
 \r
@@ -735,20 +773,25 @@ void AliFemtoK0Analysis::Exec(Option_t *)
   float qinv, q0, qx, qy, qz, qLong, qSide, qOut;              //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
 \r
   for(int i=0; i<(fEvt)->fNumV0s; i++) // Current event V0\r
   {\r
     //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("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("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
+     ((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("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("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
 \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
@@ -785,6 +828,22 @@ void AliFemtoK0Analysis::Exec(Option_t *)
         \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
+        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
+        ((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

diff --git a/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.h b/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.h
index 808887e39a5543a7b17dceb6aa56cd74067c9637..0eeed687dc4dc616f4c4feee8e6a0e0ba7cfec88 100644 (file)
--- a/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.h
+++ b/PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.h
@@ -30,7 +30,7 @@ class AliESDpid;
 class AliFemtoK0Analysis : public AliAnalysisTaskSE {
  public:
   AliFemtoK0Analysis();
-  AliFemtoK0Analysis(const char *name, bool FieldPositive = kTRUE, bool OnlineCase = kTRUE, bool MeritCase = kTRUE);
+  AliFemtoK0Analysis(const char *name, bool FieldPositive = kTRUE, bool OnlineCase = kTRUE, bool MeritCase = kTRUE, float MinDL = 0.0);
   virtual ~AliFemtoK0Analysis();
   AliFemtoK0Analysis(const AliFemtoK0Analysis&);
   AliFemtoK0Analysis& operator=(const AliFemtoK0Analysis&);
@@ -58,6 +58,7 @@ class AliFemtoK0Analysis : public AliAnalysisTaskSE {
   bool fFieldPos;
   bool fOnlineCase;
   bool fMeritCase;
+  double fMinDecayLength;
   int fEventCount;
 
   AliFemtoK0EventCollection ***fEC; //!
@@ -68,9 +69,7 @@ class AliFemtoK0Analysis : public AliAnalysisTaskSE {
   const char     *fName;
   AliAODEvent    *fAOD; //!    // AOD object
   TList          *fOutputList; //! Compact Output list
-  AliAODpidUtil  *fPidAOD; //!
-  //AliPIDResponse *fPidAOD; //!
-  AliESDpid      *fPidESD; //!
+  AliPIDResponse *fPidAOD; //!
   AliESDtrack    *fPosDaughter1;//!
   AliESDtrack    *fPosDaughter2;//!
   AliESDtrack    *fNegDaughter1;//!