X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PWG4%2FJetTasks%2FAliAnalysisHelperJetTasks.cxx;h=c6866df9250e1c28aa1b3157a233f8b61543b785;hp=20008803c67041067bb00d8b1dfecd43b03f5d4c;hb=20d01a2f89ea67653ff0136e54bff7f8a5971425;hpb=188a1ba9005eb4df24d083632c3f5b058d20aca9 diff --git a/PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx b/PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx index 20008803c67..c6866df9250 100644 --- a/PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx +++ b/PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx @@ -1,11 +1,15 @@ #include "TROOT.h" +#include "TKey.h" #include "TList.h" +#include "TSystem.h" #include "TH1F.h" #include "TProfile.h" #include "THnSparse.h" #include "TFile.h" +#include "TString.h" #include "AliMCEvent.h" +#include "AliLog.h" #include "AliAODJet.h" #include "AliStack.h" #include "AliGenEventHeader.h" @@ -14,7 +18,9 @@ #include #include #include "AliAnalysisHelperJetTasks.h" - +#include "TMatrixDSym.h" +#include "TMatrixDSymEigen.h" +#include "TVector.h" ClassImp(AliAnalysisHelperJetTasks) @@ -30,7 +36,8 @@ AliGenPythiaEventHeader* AliAnalysisHelperJetTasks::GetPythiaEventHeader(AliMCE AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast(genHeader); if (!genCocktailHeader) { - Printf("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__); + AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)"); + // AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__)); return 0; } TList* headerList = genCocktailHeader->GetHeaders(); @@ -40,7 +47,7 @@ AliGenPythiaEventHeader* AliAnalysisHelperJetTasks::GetPythiaEventHeader(AliMCE break; } if(!pythiaGenHeader){ - Printf("%s %d: PythiaHeader not found!",(char*)__FILE__,__LINE__); + AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found"); return 0; } } @@ -79,8 +86,8 @@ void AliAnalysisHelperJetTasks::PrintStack(AliMCEvent *mcEvent,Int_t iFirst,Int_ -void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &nGenJets, - AliAODJet *recJets,const Int_t &nRecJets, +void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &kGenJets, + AliAODJet *recJets,const Int_t &kRecJets, Int_t *iGenIndex,Int_t *iRecIndex, Int_t iDebug,Float_t maxDist){ @@ -95,6 +102,7 @@ void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &n // and check if there is no other rec which is closer // Caveat: Close low energy/split jets may disturb this correlation + // Idea: search in two directions generated e.g (a--e) and rec (1--3) // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet // in the end we have something like this @@ -112,17 +120,26 @@ void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &n // d c // 3 e // Only entries with "3" match from both sides + + // In case we have more jets than kmaxjets only the + // first kmaxjets are searched + // all other are -1 + // use kMaxJets for a test not to fragemnt the memory... + + for(int i = 0;i < kGenJets;++i)iGenIndex[i] = -1; + for(int j = 0;j < kRecJets;++j)iRecIndex[j] = -1; + + const int kMode = 3; - + const Int_t nGenJets = TMath::Min(kMaxJets,kGenJets); + const Int_t nRecJets = TMath::Min(kMaxJets,kRecJets); if(nRecJets==0||nGenJets==0)return; - for(int i = 0;i < nGenJets;++i)iGenIndex[i] = -1; - for(int j = 0;j < nRecJets;++j)iRecIndex[j] = -1; - - UShort_t *iFlag = new UShort_t[nGenJets*nRecJets]; + // UShort_t *iFlag = new UShort_t[nGenJets*nRecJets]; + UShort_t iFlag[kMaxJets*kMaxJets]; for(int i = 0;i < nGenJets;++i){ for(int j = 0;j < nRecJets;++j){ iFlag[i*nGenJets+j] = 0; @@ -191,9 +208,6 @@ void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &n } if(iDebug>1)printf("\n"); } - - delete [] iFlag; - } @@ -226,6 +240,7 @@ void AliAnalysisHelperJetTasks::MergeOutput(char* cFiles, char* cList){ while(in1>>cFile){ fIn[ibTotal] = TFile::Open(cFile); lIn[ibTotal] = (TList*)fIn[ibTotal]->Get(cList); + Printf("Merging file %s",cFile); if(!lIn[ibTotal]){ Printf("%s:%d No list %s found, exiting...",__FILE__,__LINE__,cList); fIn[ibTotal]->ls(); @@ -292,3 +307,309 @@ void AliAnalysisHelperJetTasks::MergeOutput(char* cFiles, char* cList){ lOut->Write(lOut->GetName(),TObject::kSingleKey); fOut->Close(); } + +Bool_t AliAnalysisHelperJetTasks::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 + // This is to called in Notify and should provide the path to the AOD/ESD file + + TString file(currFile); + fXsec = 0; + fTrials = 1; + + if(file.Contains("root_archive.zip#")){ + Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact); + Ssiz_t pos = file.Index("#",1,pos1,TString::kExact); + file.Replace(pos+1,20,""); + } + else { + // not an archive take the basename.... + file.ReplaceAll(gSystem->BaseName(file.Data()),""); + } + Printf("%s",file.Data()); + + + + + TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root + if(!fxsec){ + // next trial fetch the histgram file + fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root")); + if(!fxsec){ + // not a severe condition but inciate that we have no information + return kFALSE; + } + else{ + // find the tlist we want to be independtent of the name so use the Tkey + TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); + if(!key){ + fxsec->Close(); + return kFALSE; + } + TList *list = dynamic_cast(key->ReadObj()); + if(!list){ + fxsec->Close(); + return kFALSE; + } + fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1); + fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1); + fxsec->Close(); + } + } // no tree pyxsec.root + else { + TTree *xtree = (TTree*)fxsec->Get("Xsection"); + if(!xtree){ + fxsec->Close(); + return kFALSE; + } + UInt_t ntrials = 0; + Double_t xsection = 0; + xtree->SetBranchAddress("xsection",&xsection); + xtree->SetBranchAddress("ntrials",&ntrials); + xtree->GetEntry(0); + fTrials = ntrials; + fXsec = xsection; + fxsec->Close(); + } + return kTRUE; +} + +//___________________________________________________________________________________________________________ + +Bool_t AliAnalysisHelperJetTasks::GetEventShapes(TVector3 &n01, TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes) +{ + // *** + // Event shape calculation + // sona.pochybova@cern.ch + + const Int_t kTracks = 1000; + if(nTracks>kTracks)return kFALSE; + + //variables for thrust calculation + TVector3 pTrackPerp[kTracks]; + 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 thrust02[kTracks]; + Double_t th02 = -4; + Double_t thrust03[kTracks]; + Double_t th03 = -5; + + //Sphericity calculation variables + TMatrixDSym m(3); + Double_t s00 = 0; + Double_t s01 = 0; + Double_t s02 = 0; + + Double_t s10 = 0; + Double_t s11 = 0; + Double_t s12 = 0; + + Double_t s20 = 0; + Double_t s21 = 0; + Double_t s22 = 0; + + Double_t ptot = 0; + + Double_t c = -10; + +// +//loop for thrust calculation +// + + 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++) + { + + 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); + + } + + //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 +// + for(Int_t j = 0; j < nTracks; j++) + { + s00 = s00 + (pTrack[j].Px()*pTrack[j].Px())/pTrack[j].Mag(); + s01 = s01 + (pTrack[j].Px()*pTrack[j].Py())/pTrack[j].Mag(); + s02 = s02 + (pTrack[j].Px()*pTrack[j].Pz())/pTrack[j].Mag(); + + s10 = s10 + (pTrack[j].Py()*pTrack[j].Px())/pTrack[j].Mag(); + s11 = s11 + (pTrack[j].Py()*pTrack[j].Py())/pTrack[j].Mag(); + s12 = s12 + (pTrack[j].Py()*pTrack[j].Pz())/pTrack[j].Mag(); + + s20 = s20 + (pTrack[j].Pz()*pTrack[j].Px())/pTrack[j].Mag(); + s21 = s21 + (pTrack[j].Pz()*pTrack[j].Py())/pTrack[j].Mag(); + s22 = s22 + (pTrack[j].Pz()*pTrack[j].Pz())/pTrack[j].Mag(); + + ptot += pTrack[j].Mag(); + } + + if(ptot > 0.) + { + m(0,0) = s00/ptot; + m(0,1) = s01/ptot; + m(0,2) = s02/ptot; + + m(1,0) = s10/ptot; + m(1,1) = s11/ptot; + m(1,2) = s12/ptot; + + m(2,0) = s20/ptot; + m(2,1) = s21/ptot; + m(2,2) = s22/ptot; + + TMatrixDSymEigen eigen(m); + TVectorD eigenVal = eigen.GetEigenValues(); + + Double_t sphericity = (3/2)*(eigenVal(2)+eigenVal(1)); + eventShapes[1] = sphericity; + + Double_t aplanarity = (3/2)*(eigenVal(2)); + eventShapes[2] = aplanarity; + + c = 3*(eigenVal(0)*eigenVal(1)+eigenVal(0)*eigenVal(2)+eigenVal(1)*eigenVal(2)); + eventShapes[3] = c; + } + return kTRUE; +} + + + + //__________________________________________________________________________________________________________________________