Browse Source

Added scalar product method for V2 eta-gap method using TPC

PeterParfenov 2 years ago
parent
commit
a63081cf48
6 changed files with 192 additions and 89 deletions
  1. 30 0
      bin/fit-res.cpp
  2. 82 62
      bin/get-flow.cpp
  3. 12 9
      bin/get-recentering.cpp
  4. 13 1
      bin/get-res.cpp
  5. 45 15
      eventanalysis/mfProcessEvent.cxx
  6. 10 2
      utils/mfNamespace.h

+ 30 - 0
bin/fit-res.cpp

@@ -64,6 +64,7 @@ int main(int argc, char **argv)
 
   // Init input profiles
   TProfile *pResRecoSqrTpc[NetaGaps];
+  TProfile *pResRecoSqrSpTpc[NetaGaps];
   TProfile *pResRecoSqrMcTpc[NetaGaps];
   TProfile *pResRecoSqrFHCal;
 
@@ -73,17 +74,20 @@ int main(int argc, char **argv)
   for (int i = 0; i < NetaGaps; i++)
   {
     pResRecoSqrTpc[i] = (TProfile *)fi->Get(Form("pResRecoSqrtTpc_gap%i", i));
+    pResRecoSqrSpTpc[i] = (TProfile *)fi->Get(Form("pResRecoSqrtSpTpc_gap%i", i));
     pResRecoSqrMcTpc[i] = (TProfile *)fi->Get(Form("pResRecoSqrtMcTpc_gap%i", i));
   }
   pResRecoSqrFHCal = (TProfile *)fi->Get(Form("pResRecoSqrFHCal"));
 
   // Init output histograms
   TH1F *hResTpc[NetaGaps];
+  TH1F *hResSpTpc[NetaGaps];
   TH1F *hResMcTpc[NetaGaps];
   TH1F *hResV1FHCal;
   TH1F *hResV2FHCal;
 
   TF1 *fResTpc[NetaGaps];
+  TF1 *fResSpTpc[NetaGaps];
   TF1 *fResMcTpc[NetaGaps];
   TF1 *fResV1FHCal;
   TF1 *fResV2FHCal;
@@ -96,6 +100,12 @@ int main(int argc, char **argv)
     fResTpc[i] = new TF1(Form("fResTpc_gap%i", i),
                          "pol6",
                          0.01, 80.);
+    hResSpTpc[i] = new TH1F(Form("hResSpTpc_gap%i", i),
+                            Form("hResSpTpc_gap%i", i),
+                            NresBin, 0., 100.);
+    fResSpTpc[i] = new TF1(Form("fResSpTpc_gap%i", i),
+                         "pol6",
+                         0.01, 80.);
     hResMcTpc[i] = new TH1F(Form("hResMcTpc_gap%i", i),
                             Form("hResMcTpc_gap%i", i),
                             NresBin, 0., 100.);
@@ -130,6 +140,17 @@ int main(int argc, char **argv)
       hResTpc[iGap]->SetBinContent(iBin + 1, res);
       hResTpc[iGap]->SetBinError(iBin + 1, Eres);
     }
+    for (int iBin = 0; iBin < pResRecoSqrSpTpc[iGap]->GetNbinsX(); iBin++)
+    {
+      double res2 = pResRecoSqrSpTpc[iGap]->GetBinContent(iBin + 1);
+      double Eres2 = pResRecoSqrSpTpc[iGap]->GetBinError(iBin + 1);
+      if (res2 <= 0)
+        continue;
+      double res = Sqrt(res2);
+      double Eres = Eres2 / (2 * res2);
+      hResSpTpc[iGap]->SetBinContent(iBin + 1, res);
+      hResSpTpc[iGap]->SetBinError(iBin + 1, Eres);
+    }
 
     for (int iBin = 0; iBin < pResRecoSqrMcTpc[iGap]->GetNbinsX(); iBin++)
     {
@@ -168,6 +189,7 @@ int main(int argc, char **argv)
   for (int iGap = 0; iGap < NetaGaps; iGap++)
   {
     hResTpc[iGap]->Fit(fResTpc[iGap], "WR");
+    hResSpTpc[iGap]->Fit(fResSpTpc[iGap], "WR");
     hResMcTpc[iGap]->Fit(fResMcTpc[iGap], "WR");
   }
   hResV1FHCal->Fit(fResV1FHCal, "WR");
@@ -181,12 +203,20 @@ int main(int argc, char **argv)
   {
     hResTpc[i]->Write();
   }
+  for (int i = 0; i < NetaGaps; i++)
+  {
+    hResSpTpc[i]->Write();
+  }
   hResV1FHCal->Write();
   hResV2FHCal->Write();
   for (int i = 0; i < NetaGaps; i++)
   {
     fResTpc[i]->Write();
   }
+  for (int i = 0; i < NetaGaps; i++)
+  {
+    fResSpTpc[i]->Write();
+  }
   fResV1FHCal->Write();
   fResV2FHCal->Write();
   for (int i = 0; i < NetaGaps; i++)

+ 82 - 62
bin/get-flow.cpp

@@ -26,6 +26,8 @@ using mfEp::EpNames;
 using mfEp::NetaGaps;
 using mfEp::Nflattening;
 using mfEp::NsubEvents;
+using mfEp::SpTpcLEp;
+using mfEp::SpTpcREp;
 using mfEp::TpcEtaGap;
 using mfEp::TpcLEp;
 using mfEp::TpcREp;
@@ -160,9 +162,11 @@ int main(int argc, char **argv)
 
   // Init output profiles
   TProfile2D *pv2PtTpc[NetaGaps][Npid];
+  TProfile2D *pv2PtSpTpc[NetaGaps][Npid];
   TProfile2D *pv1PtFHCal[Npid];
   TProfile2D *pv2PtFHCal[Npid];
   TProfile2D *pv2YTpc[NetaGaps][Npid];
+  TProfile2D *pv2YSpTpc[NetaGaps][Npid];
   TProfile2D *pv1YFHCal[Npid];
   TProfile2D *pv2YFHCal[Npid];
 
@@ -191,6 +195,13 @@ int main(int argc, char **argv)
                                            Form("pv2YTpc_gap%i_pid%i", iGap, iPid),
                                            320, -1.6, 1.6, NresBin, 0., 100.);
 
+      pv2PtSpTpc[iGap][iPid] = new TProfile2D(Form("pv2PtSpTpc_gap%i_pid%i", iGap, iPid),
+                                              Form("pv2PtSpTpc_gap%i_pid%i", iGap, iPid),
+                                              300, 0., 3., NresBin, 0., 100.);
+      pv2YSpTpc[iGap][iPid] = new TProfile2D(Form("pv2YSpTpc_gap%i_pid%i", iGap, iPid),
+                                             Form("pv2YSpTpc_gap%i_pid%i", iGap, iPid),
+                                             320, -1.6, 1.6, NresBin, 0., 100.);
+
       pv2PtMcTpc[iGap][iPid] = new TProfile2D(Form("pv2PtMcTpc_gap%i_pid%i", iGap, iPid),
                                               Form("pv2PtMcTpc_gap%i_pid%i", iGap, iPid),
                                               300, 0., 3., NresBin, 0., 100.);
@@ -261,6 +272,7 @@ int main(int argc, char **argv)
   inChain->SetBranchAddress("FHCalModules", &fhcalmodules);
 
   TF1 *fResTpc[NetaGaps];
+  TF1 *fResSpTpc[NetaGaps];
   TF1 *fResMcTpc[NetaGaps];
   TF1 *fResv1FHCal, *fResv2FHCal;
 
@@ -270,6 +282,7 @@ int main(int argc, char **argv)
   for (int i = 0; i < NetaGaps; i++)
   {
     fResTpc[i] = (TF1 *)fiRes->Get(Form("fResTpc_gap%i", i));
+    fResSpTpc[i] = (TF1 *)fiRes->Get(Form("fResSpTpc_gap%i", i));
     fResMcTpc[i] = (TF1 *)fiRes->Get(Form("fResMcTpc_gap%i", i));
   }
   fResv1FHCal = (TF1 *)fiRes->Get(Form("fResV1FHCal"));
@@ -278,6 +291,8 @@ int main(int argc, char **argv)
   std::vector<float> vPsiTpcR, vPsiTpcL;
   std::vector<float> vPsiMcTpcR, vPsiMcTpcL;
   float PsiFHCal;
+  std::vector<TVector2> Q_R_SP, Q_L_SP;
+  TVector2 u_R_SP, u_L_SP;
 
   // Start event loop
   int n_entries = inChain->GetEntries();
@@ -306,6 +321,8 @@ int main(int argc, char **argv)
     mfEventInfo evInfoTpcRaw = procEvent->ProcessRecoTpcEP(recoTracks, 2.);
     vPsiTpcR.clear();
     vPsiTpcL.clear();
+    Q_R_SP.clear();
+    Q_L_SP.clear();
     PsiFHCal = 0.;
     if (evInfoTpcRaw.GetQvectors().size() == 0)
       continue;
@@ -332,6 +349,16 @@ int main(int argc, char **argv)
         vPsiTpcR.push_back(-999.);
         vPsiTpcL.push_back(-999.);
       }
+      if (selector->isGoodTpcEp(evInfoTpc.GetEpMult(SpTpcREp[iGap]), evInfoTpc.GetEpMult(SpTpcLEp[iGap])))
+      {
+        Q_R_SP.push_back((TVector2)evInfoTpc.GetQvector(SpTpcREp[iGap]));
+        Q_L_SP.push_back((TVector2)evInfoTpc.GetQvector(SpTpcLEp[iGap]));
+      }
+      else
+      {
+        Q_R_SP.push_back({0., 0.});
+        Q_L_SP.push_back({0., 0.});
+      }
     }
 
     mfEventInfo evInfoFHCalRaw = procEvent->ProcessRecoFHCalEP(fhcalmodules, 1.);
@@ -419,14 +446,28 @@ int main(int argc, char **argv)
       // v2 (TPC EP)
       for (int iGap = 0; iGap < NetaGaps; iGap++)
       {
-        double v2Tpc, v2TpcMc;
-        if (vPsiTpcR.at(iGap) == -999.)
+        double v2Tpc, v2TpcMc, v2SpTpc;
+        if (vPsiTpcR.at(iGap) == -999. || (Q_R_SP.at(iGap).X() == 0. && Q_R_SP.at(iGap).Y() == 0.))
+          continue;
+        if (vPsiTpcL.at(iGap) == -999. || (Q_L_SP.at(iGap).X() == 0. && Q_L_SP.at(iGap).Y() == 0.))
           continue;
         if (recoTrack->GetEta() < -1. * TpcEtaGap[iGap])
         {
+          u_L_SP.Set(Cos(2. * recoTrack->GetPhi()), Sin(2. * recoTrack->GetPhi()));
+          v2SpTpc = (double)(u_L_SP * Q_R_SP.at(iGap)) / fResSpTpc[iGap]->Eval(cent);
           v2Tpc = Cos(2. * (recoTrack->GetPhi() - vPsiTpcR.at(iGap))) / fResTpc[iGap]->Eval(cent);
           v2TpcMc = Cos(2. * (mcTrack->GetPhi() - vPsiTpcR.at(iGap))) / fResTpc[iGap]->Eval(cent);
+        }
+        if (recoTrack->GetEta() > TpcEtaGap[iGap])
+        {
+          u_R_SP.Set(Cos(2. * recoTrack->GetPhi()), Sin(2. * recoTrack->GetPhi()));
+          v2SpTpc = (double)(u_R_SP * Q_L_SP.at(iGap)) / fResSpTpc[iGap]->Eval(cent);
+          v2Tpc = Cos(2. * (recoTrack->GetPhi() - vPsiTpcL.at(iGap))) / fResTpc[iGap]->Eval(cent);
+          v2TpcMc = Cos(2. * (mcTrack->GetPhi() - vPsiTpcL.at(iGap))) / fResTpc[iGap]->Eval(cent);
+        }
 
+        if (Abs(recoTrack->GetEta()) > TpcEtaGap[iGap])
+        {
           if (recoTrack->GetCharge() > 0)
             pv2PtTpc[iGap][0]->Fill(recoTrack->GetPt(), cent, v2Tpc);
           if (selector->isPion(recoTrack) && recoTrack->GetCharge() > 0)
@@ -459,74 +500,21 @@ int main(int argc, char **argv)
           if (selector->isPdgProton(mcTrack)) //if (selector->isProton(recoTrack) && recoTrack->GetCharge() > 0)
             pv2YRecoMcTpc[iGap][3]->Fill(RapidityMc, cent, v2TpcMc);
 
-          if (recoTrack->GetCharge() < 0)
-            pv2PtTpc[iGap][0]->Fill(recoTrack->GetPt(), cent, v2Tpc);
-          if (selector->isPion(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2PtTpc[iGap][1]->Fill(recoTrack->GetPt(), cent, v2Tpc);
-          if (selector->isKaon(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2PtTpc[iGap][2]->Fill(recoTrack->GetPt(), cent, v2Tpc);
-          if (selector->isProton(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2PtTpc[iGap][3]->Fill(recoTrack->GetPt(), cent, v2Tpc);
-
-          if (selector->isPion(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2YTpc[iGap][1]->Fill(Rapidity, cent, v2Tpc);
-          if (selector->isKaon(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2YTpc[iGap][2]->Fill(Rapidity, cent, v2Tpc);
-          if (selector->isProton(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2YTpc[iGap][3]->Fill(Rapidity, cent, v2Tpc);
-
-          if (recoTrack->GetCharge() < 0)
-            pv2PtRecoMcTpc[iGap][4]->Fill(mcTrack->GetPt(), cent, v2TpcMc);
-          if (selector->isPdgAntiPion(mcTrack)) //if (selector->isPion(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2PtRecoMcTpc[iGap][5]->Fill(mcTrack->GetPt(), cent, v2TpcMc);
-          if (selector->isPdgAntiKaon(mcTrack)) //if (selector->isKaon(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2PtRecoMcTpc[iGap][6]->Fill(mcTrack->GetPt(), cent, v2TpcMc);
-          if (selector->isPdgAntiProton(mcTrack)) //if (selector->isProton(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2PtRecoMcTpc[iGap][7]->Fill(mcTrack->GetPt(), cent, v2TpcMc);
-
-          if (selector->isPdgAntiPion(mcTrack)) //if (selector->isPion(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2YRecoMcTpc[iGap][5]->Fill(RapidityMc, cent, v2TpcMc);
-          if (selector->isPdgAntiKaon(mcTrack)) //if (selector->isKaon(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2YRecoMcTpc[iGap][6]->Fill(RapidityMc, cent, v2TpcMc);
-          if (selector->isPdgAntiProton(mcTrack)) //if (selector->isProton(recoTrack) && recoTrack->GetCharge() < 0)
-            pv2YRecoMcTpc[iGap][7]->Fill(RapidityMc, cent, v2TpcMc);
-        }
-        if (recoTrack->GetEta() > TpcEtaGap[iGap])
-        {
-          v2Tpc = Cos(2. * (recoTrack->GetPhi() - vPsiTpcL.at(iGap))) / fResTpc[iGap]->Eval(cent);
-          v2TpcMc = Cos(2. * (mcTrack->GetPhi() - vPsiTpcL.at(iGap))) / fResTpc[iGap]->Eval(cent);
-
           if (recoTrack->GetCharge() > 0)
-            pv2PtTpc[iGap][0]->Fill(recoTrack->GetPt(), cent, v2Tpc);
+            pv2PtSpTpc[iGap][0]->Fill(recoTrack->GetPt(), cent, v2SpTpc);
           if (selector->isPion(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2PtTpc[iGap][1]->Fill(recoTrack->GetPt(), cent, v2Tpc);
+            pv2PtSpTpc[iGap][1]->Fill(recoTrack->GetPt(), cent, v2SpTpc);
           if (selector->isKaon(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2PtTpc[iGap][2]->Fill(recoTrack->GetPt(), cent, v2Tpc);
+            pv2PtSpTpc[iGap][2]->Fill(recoTrack->GetPt(), cent, v2SpTpc);
           if (selector->isProton(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2PtTpc[iGap][3]->Fill(recoTrack->GetPt(), cent, v2Tpc);
+            pv2PtSpTpc[iGap][3]->Fill(recoTrack->GetPt(), cent, v2SpTpc);
 
           if (selector->isPion(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2YTpc[iGap][1]->Fill(Rapidity, cent, v2Tpc);
+            pv2YSpTpc[iGap][1]->Fill(Rapidity, cent, v2SpTpc);
           if (selector->isKaon(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2YTpc[iGap][2]->Fill(Rapidity, cent, v2Tpc);
+            pv2YSpTpc[iGap][2]->Fill(Rapidity, cent, v2SpTpc);
           if (selector->isProton(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2YTpc[iGap][3]->Fill(Rapidity, cent, v2Tpc);
-
-          if (recoTrack->GetCharge() > 0)
-            pv2PtRecoMcTpc[iGap][0]->Fill(mcTrack->GetPt(), cent, v2TpcMc);
-          if (selector->isPdgPion(mcTrack)) //if (selector->isPion(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2PtRecoMcTpc[iGap][1]->Fill(mcTrack->GetPt(), cent, v2TpcMc);
-          if (selector->isPdgKaon(mcTrack)) //if (selector->isKaon(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2PtRecoMcTpc[iGap][2]->Fill(mcTrack->GetPt(), cent, v2TpcMc);
-          if (selector->isPdgProton(mcTrack)) //if (selector->isProton(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2PtRecoMcTpc[iGap][3]->Fill(mcTrack->GetPt(), cent, v2TpcMc);
-
-          if (selector->isPdgPion(mcTrack)) //if (selector->isPion(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2YRecoMcTpc[iGap][1]->Fill(RapidityMc, cent, v2TpcMc);
-          if (selector->isPdgKaon(mcTrack)) //if (selector->isKaon(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2YRecoMcTpc[iGap][2]->Fill(RapidityMc, cent, v2TpcMc);
-          if (selector->isPdgProton(mcTrack)) //if (selector->isProton(recoTrack) && recoTrack->GetCharge() > 0)
-            pv2YRecoMcTpc[iGap][3]->Fill(RapidityMc, cent, v2TpcMc);
+            pv2YSpTpc[iGap][3]->Fill(Rapidity, cent, v2SpTpc);
 
           if (recoTrack->GetCharge() < 0)
             pv2PtTpc[iGap][4]->Fill(recoTrack->GetPt(), cent, v2Tpc);
@@ -559,6 +547,22 @@ int main(int argc, char **argv)
             pv2YRecoMcTpc[iGap][6]->Fill(RapidityMc, cent, v2TpcMc);
           if (selector->isPdgAntiProton(mcTrack)) //if (selector->isProton(recoTrack) && recoTrack->GetCharge() < 0)
             pv2YRecoMcTpc[iGap][7]->Fill(RapidityMc, cent, v2TpcMc);
+
+          if (recoTrack->GetCharge() < 0)
+            pv2PtSpTpc[iGap][4]->Fill(recoTrack->GetPt(), cent, v2SpTpc);
+          if (selector->isPion(recoTrack) && recoTrack->GetCharge() < 0)
+            pv2PtSpTpc[iGap][5]->Fill(recoTrack->GetPt(), cent, v2SpTpc);
+          if (selector->isKaon(recoTrack) && recoTrack->GetCharge() < 0)
+            pv2PtSpTpc[iGap][6]->Fill(recoTrack->GetPt(), cent, v2SpTpc);
+          if (selector->isProton(recoTrack) && recoTrack->GetCharge() < 0)
+            pv2PtSpTpc[iGap][7]->Fill(recoTrack->GetPt(), cent, v2SpTpc);
+
+          if (selector->isPion(recoTrack) && recoTrack->GetCharge() < 0)
+            pv2YTpc[iGap][5]->Fill(Rapidity, cent, v2Tpc);
+          if (selector->isKaon(recoTrack) && recoTrack->GetCharge() < 0)
+            pv2YTpc[iGap][6]->Fill(Rapidity, cent, v2Tpc);
+          if (selector->isProton(recoTrack) && recoTrack->GetCharge() < 0)
+            pv2YTpc[iGap][7]->Fill(Rapidity, cent, v2Tpc);
         }
       }
 
@@ -905,6 +909,7 @@ int main(int argc, char **argv)
   TFile *fo = new TFile(oFileName.Data(), "recreate");
   fo->cd();
   fo->mkdir("TpcEp");
+  fo->mkdir("TpcSp");
   fo->mkdir("McTpcEp");
   fo->mkdir("RecoMcTpcEp");
   fo->mkdir("FHCalEp");
@@ -925,6 +930,21 @@ int main(int argc, char **argv)
       pv2YTpc[iGap][iPid]->Write();
     }
   }
+  fo->cd("TpcSp");
+  for (int iGap = 0; iGap < NetaGaps; iGap++)
+  {
+    for (int iPid = 0; iPid < Npid; iPid++)
+    {
+      pv2PtSpTpc[iGap][iPid]->Write();
+    }
+  }
+  for (int iGap = 0; iGap < NetaGaps; iGap++)
+  {
+    for (int iPid = 0; iPid < Npid; iPid++)
+    {
+      pv2YSpTpc[iGap][iPid]->Write();
+    }
+  }
   fo->cd("FHCalEp");
   for (int iPid = 0; iPid < Npid; iPid++)
   {

+ 12 - 9
bin/get-recentering.cpp

@@ -156,20 +156,23 @@ int main(int argc, char **argv)
     // std::cout << "\tNtracks = " << recoTracks->GetEntriesFast() << std::endl;
 
     mfEventInfo evInfoTpc = procEvent->ProcessRecoTpcEP(recoTracks, 2.);
-    if (evInfoTpc.GetQvectors().size() == 0)
-      continue;
-
-    for (int iSub = 0; iSub < NsubEventsTpc; iSub++)
+    if (evInfoTpc.GetQvectors().size() != 0)
     {
-      pQxAvg[iSub]->Fill(cent, evInfoTpc.GetQvector(iSub).X());
-      pQyAvg[iSub]->Fill(cent, evInfoTpc.GetQvector(iSub).Y());
+      for (int iSub = 0; iSub < NsubEventsTpc; iSub++)
+      {
+        pQxAvg[iSub]->Fill(cent, evInfoTpc.GetQvector(iSub).X());
+        pQyAvg[iSub]->Fill(cent, evInfoTpc.GetQvector(iSub).Y());
+      }
     }
 
     mfEventInfo evInfoFHCal = procEvent->ProcessRecoFHCalEP(fhcalmodules, 1.);
-    for (int iSub = NsubEventsTpc; iSub < NsubEvents; iSub++)
+    if (evInfoFHCal.GetQvectors().size() != 0)
     {
-      pQxAvg[iSub]->Fill(cent, evInfoFHCal.GetQvector(iSub).X());
-      pQyAvg[iSub]->Fill(cent, evInfoFHCal.GetQvector(iSub).Y());
+      for (int iSub = NsubEventsTpc; iSub < NsubEvents; iSub++)
+      {
+        pQxAvg[iSub]->Fill(cent, evInfoFHCal.GetQvector(iSub).X());
+        pQyAvg[iSub]->Fill(cent, evInfoFHCal.GetQvector(iSub).Y());
+      }
     }
   }
 

+ 13 - 1
bin/get-res.cpp

@@ -28,6 +28,8 @@ using mfEp::NsubEvents;
 using mfEp::TpcEtaGap;
 using mfEp::TpcLEp;
 using mfEp::TpcREp;
+using mfEp::SpTpcLEp;
+using mfEp::SpTpcREp;
 using mfRes::NresBin;
 
 using TMath::ATan2;
@@ -140,6 +142,7 @@ int main(int argc, char **argv)
 
   // Init output profiles
   TProfile *pResRecoSqrTpc[NetaGaps];
+  TProfile *pResRecoSqrSpTpc[NetaGaps];
   TProfile *pResRpTpc[NetaGaps];
   TProfile *pResRpTpcR[NetaGaps];
   TProfile *pResRpTpcL[NetaGaps];
@@ -160,6 +163,9 @@ int main(int argc, char **argv)
     pResRecoSqrTpc[i] = new TProfile(Form("pResRecoSqrtTpc_gap%i", i),
                                      Form("pResRecoSqrtTpc_gap%i", i),
                                      NresBin, 0., 100.);
+    pResRecoSqrSpTpc[i] = new TProfile(Form("pResRecoSqrtSpTpc_gap%i", i),
+                                       Form("pResRecoSqrtSpTpc_gap%i", i),
+                                       NresBin, 0., 100.);
     pResRpTpc[i] = new TProfile(Form("pResRpTpc_gap%i", i),
                                 Form("pResRpTpc_gap%i", i),
                                 NresBin, 0, 100);
@@ -222,6 +228,7 @@ int main(int argc, char **argv)
   mfSelector *selector = new mfSelector();
   mfInputReader *inputReader = new mfInputReader();
   inputReader->SetDcaFitFile(dcaFileName);
+  TVector2 vQ_SP_L, vQ_SP_R;
   for (int iEv = 0; iEv < n_entries; iEv++)
   {
     if (iEv % 1000 == 0)
@@ -233,7 +240,7 @@ int main(int argc, char **argv)
     if (cent == -1)
       continue;
     procEvent->SetInputReader(inputReader);
-    
+
     mfAcceptanceFilter *AccFilter = new mfAcceptanceFilter();
     if (isAcceptanceFilterOn)
       AccFilter->Activate();
@@ -257,10 +264,14 @@ int main(int argc, char **argv)
       double PsiL = evInfoTpcFlat.GetQvector(TpcLEp[iGap]).GetPhi();
 
       double resReco = Cos(2. * (PsiR - PsiL));
+      vQ_SP_L.Set(evInfoTpcFlat.GetQvector(SpTpcLEp[iGap]).X(), evInfoTpcFlat.GetQvector(SpTpcLEp[iGap]).Y());
+      vQ_SP_R.Set(evInfoTpcFlat.GetQvector(SpTpcREp[iGap]).X(), evInfoTpcFlat.GetQvector(SpTpcREp[iGap]).Y());
+      double resRecoSP = (double)(vQ_SP_R * vQ_SP_L);
       double resRpR = Cos(2. * (PsiR - mcEvent->GetPhiRP()));
       double resRpL = Cos(2. * (PsiL - mcEvent->GetPhiRP()));
 
       pResRecoSqrTpc[iGap]->Fill(cent, resReco);
+      pResRecoSqrSpTpc[iGap]->Fill(cent, resRecoSP);
       pResRpTpcR[iGap]->Fill(cent, resRpR);
       pResRpTpcL[iGap]->Fill(cent, resRpL);
       pResRpTpc[iGap]->Fill(cent, resRpR);
@@ -324,6 +335,7 @@ int main(int argc, char **argv)
   for (int iGap = 0; iGap < NetaGaps; iGap++)
   {
     pResRecoSqrTpc[iGap]->Write();
+    pResRecoSqrSpTpc[iGap]->Write();
     pResRpTpcR[iGap]->Write();
     pResRpTpcL[iGap]->Write();
     pResRpTpc[iGap]->Write();

+ 45 - 15
eventanalysis/mfProcessEvent.cxx

@@ -17,6 +17,8 @@ using mfEp::NetaGaps;
 using mfEp::Nmodules;
 using mfEp::NsubEvents;
 using mfEp::NsubEventsTpc;
+using mfEp::SpTpcLEp;
+using mfEp::SpTpcREp;
 using mfEp::TpcEtaGap;
 using mfEp::TpcLEp;
 using mfEp::TpcREp;
@@ -106,7 +108,11 @@ mfEventInfo mfProcessEvent::ProcessRecoTpcEP(TClonesArray *const &recoTracks, Do
         Qx.at(TpcLEp[iGap]) += recoTrack->GetPt() * Cos(harmonic * recoTrack->GetPhi());
         Qy.at(TpcLEp[iGap]) += recoTrack->GetPt() * Sin(harmonic * recoTrack->GetPhi());
         Wt.at(TpcLEp[iGap]) += recoTrack->GetPt();
+        Qx.at(SpTpcLEp[iGap]) += 1. * Cos(harmonic * recoTrack->GetPhi());
+        Qy.at(SpTpcLEp[iGap]) += 1. * Sin(harmonic * recoTrack->GetPhi());
+        Wt.at(SpTpcLEp[iGap]) += 1.;
         EpMult.at(TpcLEp[iGap]) += 1;
+        EpMult.at(SpTpcLEp[iGap]) += 1;
       }
       // Tpc right
       if (recoTrack->GetEta() > TpcEtaGap[iGap])
@@ -114,15 +120,18 @@ mfEventInfo mfProcessEvent::ProcessRecoTpcEP(TClonesArray *const &recoTracks, Do
         Qx.at(TpcREp[iGap]) += recoTrack->GetPt() * Cos(harmonic * recoTrack->GetPhi());
         Qy.at(TpcREp[iGap]) += recoTrack->GetPt() * Sin(harmonic * recoTrack->GetPhi());
         Wt.at(TpcREp[iGap]) += recoTrack->GetPt();
+        Qx.at(SpTpcREp[iGap]) += 1. * Cos(harmonic * recoTrack->GetPhi());
+        Qy.at(SpTpcREp[iGap]) += 1. * Sin(harmonic * recoTrack->GetPhi());
+        Wt.at(SpTpcREp[iGap]) += 1.;
         EpMult.at(TpcREp[iGap]) += 1;
+        EpMult.at(SpTpcREp[iGap]) += 1;
       }
     }
   }
-  for (int i = 0; i < NsubEventsTpc; i++)
+  for (int i = EpNames::kTpc_R_EtaGap1; i <= EpNames::kTpc_L_EtaGap3; i++)
   {
     if (Wt.at(i) == 0 || EpMult.at(i) == 0)
     {
-      // std::cerr << "No tracks for EP" << std::endl;
       return fEventInfo;
     }
     else
@@ -130,8 +139,17 @@ mfEventInfo mfProcessEvent::ProcessRecoTpcEP(TClonesArray *const &recoTracks, Do
       Qx.at(i) /= Wt.at(i);
       Qy.at(i) /= Wt.at(i);
       Qvectors.at(i).Set(Qx.at(i), Qy.at(i));
-      // std::cout << "Qx[" << i << "] = " << Qx.at(i) << ", "
-      //           << "Qy[" << i << "] = " << Qy.at(i) << std::endl;
+    }
+  }
+  for (int i = EpNames::kSpTpc_R_EtaGap1; i <= EpNames::kSpTpc_L_EtaGap3; i++)
+  {
+    if (Wt.at(i) == 0 || EpMult.at(i) == 0)
+    {
+      return fEventInfo;
+    }
+    else
+    {
+      Qvectors.at(i).Set(Qx.at(i), Qy.at(i));
     }
   }
   fEventInfo.SetMultiplicity(mult);
@@ -178,10 +196,6 @@ mfEventInfo mfProcessEvent::ProcessRecoFHCalEP(TClonesArray *const &fhcalModules
 
     double phi = GetFHCalPhi(iMod);
 
-    // std::cout << "\tModule " << iMod
-    //           << ": phi " << phi
-    //           << ", Energy " << fhcalModule->GetEnergy() << std::endl;
-
     int w_sign = (iMod < Nmodules) ? 1 : -1;
     // for left FHCal modules
     if (iMod < Nmodules)
@@ -203,7 +217,7 @@ mfEventInfo mfProcessEvent::ProcessRecoFHCalEP(TClonesArray *const &fhcalModules
     Wt.at(EpNames::kFHCal) += fhcalModule->GetEnergy();
   }
 
-  for (int i = NsubEventsTpc; i < NsubEvents; i++)
+  for (int i = EpNames::kFHCal_R; i <= EpNames::kFHCal; i++)
   {
     if (Wt.at(i) == 0)
     {
@@ -215,8 +229,6 @@ mfEventInfo mfProcessEvent::ProcessRecoFHCalEP(TClonesArray *const &fhcalModules
       Qx.at(i) /= Wt.at(i);
       Qy.at(i) /= Wt.at(i);
       Qvectors.at(i).Set(Qx.at(i), Qy.at(i));
-      // std::cout << "Qx[" << i << "] = " << Qx.at(i) << ", "
-      //           << "Qy[" << i << "] = " << Qy.at(i) << std::endl;
     }
   }
 
@@ -280,7 +292,12 @@ mfEventInfo mfProcessEvent::ProcessMcTpcEP(TClonesArray *const &mcTracks, Double
         Qx.at(TpcLEp[iGap]) += mcTrack->GetPt() * Cos(harmonic * mcTrack->GetPhi());
         Qy.at(TpcLEp[iGap]) += mcTrack->GetPt() * Sin(harmonic * mcTrack->GetPhi());
         Wt.at(TpcLEp[iGap]) += mcTrack->GetPt();
+        Qx.at(SpTpcLEp[iGap]) += 1. * Cos(harmonic * mcTrack->GetPhi());
+        Qy.at(SpTpcLEp[iGap]) += 1. * Sin(harmonic * mcTrack->GetPhi());
+        Wt.at(TpcLEp[iGap]) += 1.;
+        Wt.at(SpTpcLEp[iGap]) += 1.;
         EpMult.at(TpcLEp[iGap]) += 1;
+        EpMult.at(SpTpcLEp[iGap]) += 1;
       }
       // Tpc right
       if (mcTrack->GetEta() > TpcEtaGap[iGap])
@@ -288,15 +305,19 @@ mfEventInfo mfProcessEvent::ProcessMcTpcEP(TClonesArray *const &mcTracks, Double
         Qx.at(TpcREp[iGap]) += mcTrack->GetPt() * Cos(harmonic * mcTrack->GetPhi());
         Qy.at(TpcREp[iGap]) += mcTrack->GetPt() * Sin(harmonic * mcTrack->GetPhi());
         Wt.at(TpcREp[iGap]) += mcTrack->GetPt();
+        Qx.at(SpTpcREp[iGap]) += 1. * Cos(harmonic * mcTrack->GetPhi());
+        Qy.at(SpTpcREp[iGap]) += 1. * Sin(harmonic * mcTrack->GetPhi());
+        Wt.at(SpTpcREp[iGap]) += 1.;
+        Wt.at(SpTpcREp[iGap]) += 1.;
         EpMult.at(TpcREp[iGap]) += 1;
+        EpMult.at(SpTpcREp[iGap]) += 1;
       }
     }
   }
-  for (int i = 0; i < NsubEventsTpc; i++)
+  for (int i = EpNames::kTpc_R_EtaGap1; i <= EpNames::kTpc_L_EtaGap3; i++)
   {
     if (Wt.at(i) == 0 || EpMult.at(i) == 0)
     {
-      // std::cerr << "No tracks for EP" << std::endl;
       return fEventInfo;
     }
     else
@@ -304,8 +325,17 @@ mfEventInfo mfProcessEvent::ProcessMcTpcEP(TClonesArray *const &mcTracks, Double
       Qx.at(i) /= Wt.at(i);
       Qy.at(i) /= Wt.at(i);
       Qvectors.at(i).Set(Qx.at(i), Qy.at(i));
-      // std::cout << "Qx[" << i << "] = " << Qx.at(i) << ", "
-      //           << "Qy[" << i << "] = " << Qy.at(i) << std::endl;
+    }
+  }
+  for (int i = EpNames::kSpTpc_R_EtaGap1; i <= EpNames::kSpTpc_L_EtaGap3; i++)
+  {
+    if (Wt.at(i) == 0 || EpMult.at(i) == 0)
+    {
+      return fEventInfo;
+    }
+    else
+    {
+      Qvectors.at(i).Set(Qx.at(i), Qy.at(i));
     }
   }
   // fEventInfo.SetMultiplicity(mult);

+ 10 - 2
utils/mfNamespace.h

@@ -42,15 +42,23 @@ namespace mfEp
     kTpc_L_EtaGap1,
     kTpc_L_EtaGap2,
     kTpc_L_EtaGap3,
+    kSpTpc_R_EtaGap1,
+    kSpTpc_R_EtaGap2,
+    kSpTpc_R_EtaGap3,
+    kSpTpc_L_EtaGap1,
+    kSpTpc_L_EtaGap2,
+    kSpTpc_L_EtaGap3,
     kFHCal_R,
     kFHCal_L,
-    kFHCal
+    kFHCal,
   };
   const int NsubEvents = kFHCal + 1;
-  const int NsubEventsTpc = kTpc_L_EtaGap3 + 1;
+  const int NsubEventsTpc = kSpTpc_L_EtaGap3 + 1;
   const float TpcEtaGap[NetaGaps] = {0.05, 0.1, 0.25};
   const int TpcREp[NetaGaps] = {kTpc_R_EtaGap1, kTpc_R_EtaGap2, kTpc_R_EtaGap3};
   const int TpcLEp[NetaGaps] = {kTpc_L_EtaGap1, kTpc_L_EtaGap2, kTpc_L_EtaGap3};
+  const int SpTpcREp[NetaGaps] = {kSpTpc_R_EtaGap1, kSpTpc_R_EtaGap2, kSpTpc_R_EtaGap3};
+  const int SpTpcLEp[NetaGaps] = {kSpTpc_L_EtaGap1, kSpTpc_L_EtaGap2, kSpTpc_L_EtaGap3};
 
   const int Nmodules = 45;
   double GetFHCalPhi(int iModule);