account for limits of thrust variable (Sona)

diff --git a/PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx b/PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx
index 6dea29c3856fc6f9b4ce3b77d570886c12e7df25..c6866df9250e1c28aa1b3157a233f8b61543b785 100644 (file)
--- a/PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx
+++ b/PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx
@@ -374,54 +374,6 @@ Bool_t AliAnalysisHelperJetTasks::PythiaInfoFromFile(const char* currFile,Float_
   return kTRUE;
 }
 
-
-Bool_t AliAnalysisHelperJetTasks::GetThrustAxis(TVector3 &n01, TVector3 * pTrack, const Int_t &nTracks)
-{
-  //
-  // fetch the thrust axis
-  // sona.pochybova@cern.ch
-
-  const Int_t kTracks = 1000;
-  if(nTracks>kTracks)return kFALSE;
-
-  TVector3 psum;
-  Double_t psum1 = 0;
-  Double_t psum2 = 0;
-  Double_t thrust[kTracks];
-  Double_t th = -3;
-  Double_t tpom = -1;
-  Int_t j = 0;
-
-  
-  //  for(Int_t j = 0; j < nTracks; j++)
-  while(TMath::Abs(th-tpom) > 10e-7 && j < nTracks)
-    {
-      th = tpom;
-      psum.SetXYZ(0., 0., 0.);
-      psum1 = 0;
-      psum2 = 0;
-      for(Int_t i = 0; i < nTracks; i++)
-       {
-         psum1 += (TMath::Abs(n01.Dot(pTrack[i])));
-         psum2 += pTrack[i].Mag();
-         
-         if (n01.Dot(pTrack[i]) > 0) psum += pTrack[i];
-         if (n01.Dot(pTrack[i]) < 0) psum -= pTrack[i];
-       }
-      
-      thrust[j] = psum1/psum2;
-      
-      //      if(th == thrust[j]) 
-      //       break;
-      
-      tpom = thrust[j];
-      
-      n01 = psum.Unit();
-      j++;
-    }
-  return kTRUE;
-}
-
 //___________________________________________________________________________________________________________
 
 Bool_t AliAnalysisHelperJetTasks::GetEventShapes(TVector3 &n01, TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes)
@@ -433,15 +385,26 @@ Bool_t AliAnalysisHelperJetTasks::GetEventShapes(TVector3 &n01, TVector3 * pTrac
   const Int_t kTracks = 1000;
   if(nTracks>kTracks)return kFALSE;
 
-  TVector3 psum;
+  //variables for thrust calculation
   TVector3 pTrackPerp[kTracks];
-  Double_t psum1 = 0;
   Double_t psum2 = 0;
+
+  TVector3 psum;
+  TVector3 psum02;
+  TVector3 psum03;
+
+  Double_t psum1 = 0;
+  Double_t psum102 = 0;
+  Double_t psum103 = 0;
+
   Double_t thrust[kTracks];
   Double_t th = -3;
-  Double_t tpom = -1;
+  Double_t thrust02[kTracks];
+  Double_t th02 = -4;
+  Double_t thrust03[kTracks];
+  Double_t th03 = -5;
 
-  //Sphericity calculation
+  //Sphericity calculation variables
   TMatrixDSym m(3);
   Double_t s00 = 0;
   Double_t s01 = 0;
@@ -462,33 +425,141 @@ Bool_t AliAnalysisHelperJetTasks::GetEventShapes(TVector3 &n01, TVector3 * pTrac
 //
 //loop for thrust calculation  
 //
-  Int_t k = 0;
-  while(TMath::Abs(th-tpom) > 10e-7 && k < nTracks)
+    
+  for(Int_t i = 0; i < nTracks; i++)
+    {
+      pTrackPerp[i].SetXYZ(pTrack[i].X(), pTrack[i].Y(), 0);
+      psum2 += pTrackPerp[i].Mag();
+    }
+
+  //additional starting axis    
+  TVector3 n02;
+  n02 = pTrack[1].Unit();
+  n02.SetZ(0.);   
+  TVector3 n03;
+  n03 = pTrack[2].Unit();
+  n03.SetZ(0.);
+
+  //switches for calculating thrust for different starting points
+  Int_t switch1 = 1;
+  Int_t switch2 = 1;
+  Int_t switch3 = 1;
+
+  //indexes for iteration of different starting points
+  Int_t l1 = 0;
+  Int_t l2 = 0;
+  Int_t l3 = 0;
+
+  //maximal number of iterations
+  //  Int_t nMaxIter = 100;
+ 
+  for(Int_t k = 0; k < nTracks; k++)
     {  
-      th = tpom;
-      psum.SetXYZ(0., 0., 0.);
-      psum1 = 0;
-      psum2 = 0;
-      for(Int_t i = 0; i < nTracks; i++)
-       {
-         pTrackPerp[i].SetXYZ(pTrack[i].X(), pTrack[i].Y(), 0);
-
-         psum1 += (TMath::Abs(n01.Dot(pTrackPerp[i])));
-         psum2 += pTrackPerp[i].Mag();
-
-         if (n01.Dot(pTrackPerp[i]) > 0) psum += pTrackPerp[i];
-         if (n01.Dot(pTrackPerp[i]) < 0) psum -= pTrackPerp[i];
-       }
       
-      thrust[k] = psum1/psum2;
+      if(switch1 == 1){
+       psum.SetXYZ(0., 0., 0.);
+       psum1 = 0;
+       for(Int_t i = 0; i < nTracks; i++)
+         {
+           psum1 += (TMath::Abs(n01.Dot(pTrackPerp[i])));
+           if (n01.Dot(pTrackPerp[i]) > 0) psum += pTrackPerp[i];
+           if (n01.Dot(pTrackPerp[i]) < 0) psum -= pTrackPerp[i];
+         }
+       thrust[l1] = psum1/psum2;
+      }
+
+      if(switch2 == 1){
+       psum02.SetXYZ(0., 0., 0.);
+       psum102 = 0;
+       for(Int_t i = 0; i < nTracks; i++)
+         {
+           psum102 += (TMath::Abs(n02.Dot(pTrackPerp[i])));
+           if (n02.Dot(pTrackPerp[i]) > 0) psum02 += pTrackPerp[i];
+           if (n02.Dot(pTrackPerp[i]) < 0) psum02 -= pTrackPerp[i];
+         }
+       thrust02[l2] = psum102/psum2;
+      }
+      
+      if(switch3 == 1){
+       psum03.SetXYZ(0., 0., 0.);
+       psum103 = 0;
+       for(Int_t i = 0; i < nTracks; i++)
+         {
+           psum103 += (TMath::Abs(n03.Dot(pTrackPerp[i])));
+           if (n03.Dot(pTrackPerp[i]) > 0) psum03 += pTrackPerp[i];
+           if (n03.Dot(pTrackPerp[i]) < 0) psum03 -= pTrackPerp[i];
+         }
+       thrust03[l3] = psum103/psum2;
+      }
+
+      //check whether thrust value converged    
+      if(TMath::Abs(th-thrust[l1]) < 10e-7){
+       switch1 = 0;
+      }
+      
+      if(TMath::Abs(th02-thrust02[l2]) < 10e-7){
+       switch2 = 0;
+      }
+
+      if(TMath::Abs(th03-thrust03[l3]) < 10e-7){
+       switch3 = 0;
+      }
+
+      //if it didn't, continue with the calculation
+      if(switch1 == 1){
+       th = thrust[l1]; 
+       n01 = psum.Unit();
+       l1++;
+      }
+
+      if(switch2 == 1){
+       th02 = thrust02[l2];
+       n02 = psum02.Unit();
+       l2++;
+      }
+
+      if(switch3 == 1){
+       th03 = thrust03[l3];
+       n03 = psum03.Unit();
+       l3++;
+      }
+
+      //if thrust values for all starting direction converged check if to the same value
+      if(switch2 == 0 && switch1 == 0 && switch3 == 0){
+       if(TMath::Abs(th-th02) < 10e-7 && TMath::Abs(th-th03) < 10e-7 && TMath::Abs(th02-th03) < 10e-7){
+         eventShapes[0] = th;
+         Printf("===== THRUST VALUE FOUND AT %d :: %f\n", k, th);
+         break;
+       }
+       //if they did not, reset switches
+       else{
+         switch1 = 1;
+         //      th = -1.;
+         switch2 = 1;
+         //      th02 = -2.;
+         switch3 = 1;
+         //      th03 = -4.;
+       }
+      }
+
+      //      Printf("========== %d +++ th :: %f=============\n", l1, th);
+      //      Printf("========== %d +++ th2 :: %f=============\n", l2, th02);
+      //      Printf("========== %d +++ th3 :: %f=============\n", l3, th03);
       
-      tpom = thrust[k];
-      n01 = psum.Unit();
-      k++;
-//      Printf("========== %d +++ th :: %f, tpom :: %f=============\n", k, th, tpom);
     }
 
-  eventShapes[0] = th;
+  //if no common limitng value was found, take the maximum and take the corresponding thrust axis
+  if(switch1 == 1 && switch2 == 1 && switch3 == 1){
+    eventShapes[0] = TMath::Max(thrust[l1-1], thrust02[l2-1]);
+    eventShapes[0] = TMath::Max(eventShapes[0], thrust03[l3-1]);
+    if(TMath::Abs(eventShapes[0]-thrust[l1-1]) < 10e-7)
+      n01 = n01;
+    if(TMath::Abs(eventShapes[0]-thrust02[l2-1]) < 10e-7)
+      n01 = n02;
+    if(TMath::Abs(eventShapes[0]-thrust03[l3-1]) < 10e-7)
+      n01 = n03;
+    Printf("NO LIMITING VALUE FOUND :: MAXIMUM = %f\n", eventShapes[0]);
+  }
 
 //
 //other event shapes variables

diff --git a/PWG4/JetTasks/AliAnalysisHelperJetTasks.h b/PWG4/JetTasks/AliAnalysisHelperJetTasks.h
index e73de686ae601b55af76836ddee0419497dbe95d..e1033e818ce75d4e31f287421142a22d04046844 100644 (file)
--- a/PWG4/JetTasks/AliAnalysisHelperJetTasks.h
+++ b/PWG4/JetTasks/AliAnalysisHelperJetTasks.h
@@ -27,7 +27,6 @@ class AliAnalysisHelperJetTasks : public TObject {
 
   static void MergeOutput(char* cFiles, char* cList = "pwg4spec"); // Merges the files in the input text file  needs the two histograms fh1PtHard_Trials, fh1Xsec and the name of the input list
   static Bool_t PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials);// get the cross section and the trails either from pyxsec.root or from pysec_hists.root
-  static Bool_t GetThrustAxis(TVector3 &n01, TVector3 * pTrack, const Int_t &nTracks);
   static Bool_t GetEventShapes(TVector3 &n01, TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes);
   enum {kMaxJets = 6}; //  needed for array size not to fragemnt memory on the heap by many new/delete 
   private:

diff --git a/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx b/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx
index 58fe1e516db15465a75ec5dfc725491a6395db77..6c44c85ce844759955f42837cab490cd1a59e3a5 100644 (file)
--- a/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx
+++ b/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx
@@ -21,7 +21,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-
+#include <cstdlib>
+#include <iostream>
 #include <TFile.h>
 #include <TH1F.h>
 #include <TH2F.h>
@@ -124,6 +125,7 @@ ClassImp(AliAnalysisTaskThreeJets)
 
                                                         fhdPhiThrustRec(0x0),
                                                         fhdPhiThrustRecALL(0x0)
+
 {
 
 }
@@ -726,23 +728,40 @@ void AliAnalysisTaskThreeJets::UserExec(Option_t * )
       pTav[i] = pTsum/nInJet[i]; 
     }
   
-  TMath::Sort(nAccTr, pTracksMC, idxTracksMC);
-  for(Int_t i = 0; i < nAccTr; i++)
+  TMath::Sort(nAllTracksMC, pTrackMCAll, idxTracksMC);
+  for(Int_t i = 0; i < nAllTracksMC; i++)
     {
-      jetTracksSortMC[i] = jetTracksMC[idxTracksMC[i]];
+      jetTracksSortMC[i] = vTrackMCAll[idxTracksMC[i]];
       pTrackMC[i].SetXYZ(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz());
       vTrackMC[i].SetPxPyPzE(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz(), jetTracksSortMC[i].E());
     }
 
   TVector3 n01MC = pTrackMC[0].Unit();
+  n01MC.SetZ(0.);
+
   //Thrust calculation, iterative method
   if(nGenSel > 1)
     { 
 //       if(fGlobVar == 1)
 //     {
-      AliAnalysisHelperJetTasks::GetEventShapes(n01MC, pTrackMC, nAccTr, eventShapes);
+      Printf("**************Shapes for MC*************");
+      AliAnalysisHelperJetTasks::GetEventShapes(n01MC, pTrackMC, nAllTracksMC, eventShapes);
 //     }
-      
+      if(eventShapes[0] < 2/TMath::Pi()){
+       Double_t eventShapesTest[4];
+       TVector3 n01Test;
+       Int_t rnd_max = nAllTracksMC;
+       Int_t k = (rand()%rnd_max)+3;
+       while(TMath::Abs(pTrackMC[k].X()) < 10e-5 && TMath::Abs(pTrackMC[k].Y()) < 10e-5){
+         k--;
+       }
+       n01Test = pTrackMC[k].Unit();
+       n01Test.SetZ(0.);
+       AliAnalysisHelperJetTasks::GetEventShapes(n01Test, pTrackMC, nAllTracksMC, eventShapesTest);
+       eventShapes[0] = TMath::Max(eventShapes[0], eventShapesTest[0]);
+       if(TMath::Abs(eventShapes[0]-eventShapesTest[0]) < 10e-7) n01MC = n01Test;
+      }
+
       Double_t s = eventShapes[1];
       Double_t a = eventShapes[2];
       Double_t c = eventShapes[3];
@@ -785,13 +804,13 @@ void AliAnalysisTaskThreeJets::UserExec(Option_t * )
       pGen[i].SetXYZ(vGen[i].Px(), vGen[i].Py(), vGen[i].Pz());
       vsumGen += vGen[i];
     }
-  if(eventShapes[0] >0.8 )
+  if(eventShapes[0] > 0.8 && nGenSel > 1)
     {
       for(Int_t i = 0; i < nGenSel; i++)
        fhdPhiThrustGen->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
     }
  
-  if(eventShapes[0] <= 0.8)   
+  if(eventShapes[0] <= 0.8 && nGenSel > 1)   
     {
       for(Int_t i = 0; i < nGenSel; i++)
        fhdPhiThrustGenALL->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
@@ -903,9 +922,9 @@ void AliAnalysisTaskThreeJets::UserExec(Option_t * )
       eRec[rj] = recSelJets[rj].E();
     }
  
-  //  Int_t nAODtracks = fAOD->GetNumberOfTracks();
+  Int_t nAODtracks = fAOD->GetNumberOfTracks();
   Int_t nTracks = 0; //tracks accepted in the whole event
-  //  Int_t nTracksALL = 0;
+  Int_t nTracksALL = 0;
   TLorentzVector jetTracks[kTracks];
   TLorentzVector jetTracksSort[kTracks];
   Double_t * eTracks = new Double_t[kTracks];
@@ -935,10 +954,10 @@ void AliAnalysisTaskThreeJets::UserExec(Option_t * )
          AliAODTrack * track = (AliAODTrack*)jetTracksAOD->At(i);
          track->GetCovarianceXYZPxPyPz(cv);
          if(cv[14] > 1000.) continue;
-         jetTrack[nTracks] = *track;
-         jetTracks[nTracks].SetPxPyPzE(jetTrack[nTracks].Px(), jetTrack[nTracks].Py(), jetTrack[nTracks].Pz(), jetTrack[nTracks].E());
-         eTracks[nTracks] = jetTracks[nTracks].E();
-         pTracks[nTracks] = jetTracks[nTracks].Pt();
+         //      jetTrack[nTracks] = *track;
+         //      jetTracks[nTracks].SetPxPyPzE(jetTrack[nTracks].Px(), jetTrack[nTracks].Py(), jetTrack[nTracks].Pz(), jetTrack[nTracks].E());
+         //      eTracks[nTracks] = jetTracks[nTracks].E();
+         //      pTracks[nTracks] = jetTracks[nTracks].Pt();
          nTracks++;
          nJetTracks++;
        }
@@ -947,12 +966,26 @@ void AliAnalysisTaskThreeJets::UserExec(Option_t * )
       jetMult[nAccJets] = jetTracksAOD->GetEntries();
       nAccJets++;
     }
-  
+ 
   if (nAccJets == 0) return;
 
+  for(Int_t i = 0; i < nAODtracks; i++)
+    {
+      AliAODTrack * track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(i));
+      if(!track) continue;
+      track->GetCovarianceXYZPxPyPz(cv);
+      if(cv[14] > 1000.) continue;
+      jetTrack[nTracksALL] = *track;
+      jetTracks[nTracksALL].SetPxPyPzE(jetTrack[nTracksALL].Px(), jetTrack[nTracksALL].Py(), jetTrack[nTracksALL].Pz(), jetTrack[nTracksALL].E());
+      eTracks[nTracksALL] = jetTracks[nTracksALL].E();
+      pTracks[nTracksALL] = jetTracks[nTracksALL].Pt();
+      nTracksALL++;
+    }
+      
+
   TLorentzVector vTrack[kTracks];
-  TMath::Sort(nTracks, pTracks, idxTracks);
-  for(Int_t i = 0; i < nTracks; i++)
+  TMath::Sort(nTracksALL, pTracks, idxTracks);
+  for(Int_t i = 0; i < nTracksALL; i++)
     {
       jetTracksSort[i] = jetTracks[idxTracks[i]];
       pTrack[i].SetXYZ(jetTracksSort[i].Px(), jetTracksSort[i].Py(), jetTracksSort[i].Pz());
@@ -968,16 +1001,35 @@ void AliAnalysisTaskThreeJets::UserExec(Option_t * )
 
   //Thrust, iterative method, AODs
   TVector3 n01 = pTrack[0].Unit();
+  n01.SetZ(0.);
   if(nAccJets > 1)
     {
 //       if(fGlobVar == 1)
 //     {
-      AliAnalysisHelperJetTasks::GetEventShapes(n01, pTrack, nTracks, eventShapesRec);
-//     }
-//       fGlobVar = 0;      
-//       Double_t Max3 = TMath::Max(eventShapesRec0[0],eventShapesRec1[0]);
-//       Double_t Max4 = TMath::Max(eventShapesRec3[0],eventShapesRec2[0]);
- 
+      Printf("*********Shapes for AOD********");
+      AliAnalysisHelperJetTasks::GetEventShapes(n01, pTrack, nTracksALL, eventShapesRec);
+      //       }
+      if(eventShapesRec[0] < 2/TMath::Pi()){
+       Double_t eventShapesTest[4];
+       TVector3 n01Test;
+       Int_t rnd_max = nTracksALL;
+       Int_t k = (rand()%rnd_max)+3;
+       while(TMath::Abs(pTrack[k].X()) < 10e-5 && TMath::Abs(pTrack[k].Y()) < 10e-5){
+         k--;
+       }
+       
+       n01Test = pTrack[k].Unit();
+       n01Test.SetZ(0.);
+       AliAnalysisHelperJetTasks::GetEventShapes(n01Test, pTrack, nTracksALL, eventShapesTest);
+       eventShapesRec[0] = TMath::Max(eventShapesRec[0], eventShapesTest[0]);
+       if(TMath::Abs(eventShapesRec[0]-eventShapesTest[0]) < 10e-7) n01 = n01Test;
+      }
+    
+       
+      //       fGlobVar = 0;      
+      //       Double_t Max3 = TMath::Max(eventShapesRec0[0],eventShapesRec1[0]);
+      //       Double_t Max4 = TMath::Max(eventShapesRec3[0],eventShapesRec2[0]);
+      
       Double_t thrust = eventShapesRec[0];
       
       if(eventShapesRec[0] > 0.8)
@@ -1020,7 +1072,6 @@ void AliAnalysisTaskThreeJets::UserExec(Option_t * )
          }
          break;
        }
-      
     }
 
   //rest frame for reconstructed jets
@@ -1171,7 +1222,7 @@ Bool_t AliAnalysisTaskThreeJets::IsPrimChar(TParticle* aParticle, Int_t aTotalPr
   {
     if (adebug)
       printf("Dropping particle because it is not charged.\n");
-    return kFALSE;
+    return kTRUE;
   }
 
   return kTRUE;