From 20d01a2f89ea67653ff0136e54bff7f8a5971425 Mon Sep 17 00:00:00 2001 From: kleinb Date: Thu, 29 Oct 2009 08:55:03 +0000 Subject: [PATCH] account for limits of thrust variable (Sona) --- PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx | 219 +++++++++++++------- PWG4/JetTasks/AliAnalysisHelperJetTasks.h | 1 - PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx | 101 ++++++--- 3 files changed, 221 insertions(+), 100 deletions(-) diff --git a/PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx b/PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx index 6dea29c3856..c6866df9250 100644 --- 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 e73de686ae6..e1033e818ce 100644 --- 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 58fe1e516db..6c44c85ce84 100644 --- a/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx +++ b/PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx @@ -21,7 +21,8 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ - +#include +#include #include #include #include @@ -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(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; -- 2.18.4