added merging of different pt_hard bins in macro and Notify() in task to fill xsection
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 15 Jan 2009 14:04:05 +0000 (14:04 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 15 Jan 2009 14:04:05 +0000 (14:04 +0000)
PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx
PWG4/JetTasks/AliAnalysisTaskJetSpectrum.h
PWG4/macros/runJetSpectrumUnfolding.C

index 26948e55ae21b654e755a80a1d326a365c7c5fe9..e302e7e63d8b8bc681c62f8517526ba2923697ff 100644 (file)
@@ -23,6 +23,7 @@
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TH3F.h>
+#include <TProfile.h>
 #include <TList.h>
 #include <TLorentzVector.h>
 #include <TClonesArray.h>
@@ -68,6 +69,8 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
   fLimitGenJetEta(kFALSE),
   fAnalysisType(0),
   fExternalWeight(1),    
+  fh1Xsec(0x0),
+  fh1Trials(0x0),
   fh1PtHard(0x0),
   fh1PtHard_NoW(0x0),  
   fh1PtHard_Trials(0x0),
@@ -108,6 +111,8 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
   fLimitGenJetEta(kFALSE),
   fAnalysisType(0),
   fExternalWeight(1),    
+  fh1Xsec(0x0),
+  fh1Trials(0x0),
   fh1PtHard(0x0),
   fh1PtHard_NoW(0x0),  
   fh1PtHard_Trials(0x0),
@@ -139,6 +144,52 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
 
 
 
+Bool_t AliAnalysisTaskJetSpectrum::Notify()
+{
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  Double_t xsection = 0;
+  UInt_t   ntrials  = 0;
+  if(tree){
+    TFile *curfile = tree->GetCurrentFile();
+    if (!curfile) {
+      Error("Notify","No current file");
+      return kFALSE;
+    }
+    if(!fh1Xsec||!fh1Trials){
+      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
+      return kFALSE;
+    }
+
+    TString fileName(curfile->GetName());
+    if(fileName.Contains("AliESDs.root")){
+        fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
+    }
+    else if(fileName.Contains("AliAOD.root")){
+        fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
+    }
+    else if(fileName.Contains("galice.root")){
+        // for running with galice and kinematics alone...                      
+        fileName.ReplaceAll("galice.root", "pyxsec.root");
+    }
+    TFile *fxsec = TFile::Open(fileName.Data());
+    if(!fxsec){
+      Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
+      // no a severe condition
+      return kTRUE;
+    }
+    TTree *xtree = (TTree*)fxsec->Get("Xsection");
+    if(!xtree){
+      Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+    }
+    xtree->SetBranchAddress("xsection",&xsection);
+    xtree->SetBranchAddress("ntrials",&ntrials);
+    xtree->GetEntry(0);
+    fh1Xsec->Fill("<#sigma>",xsection);
+    fh1Trials->Fill("#sum{ntrials}",ntrials);
+  }
+  
+  return kTRUE;
+}
 
 void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
 {
@@ -226,6 +277,12 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
   const Int_t nBinFrag = 25;
 
 
+  fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+
+  fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
+  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+
   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 
   fh1PtHard_NoW = new TH1F("fh1PtHard_NoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
@@ -307,6 +364,8 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
   }
   const Int_t saveLevel = 1; // large save level more histos
   if(saveLevel>0){
+    fHistList->Add(fh1Xsec);
+    fHistList->Add(fh1Trials);
     fHistList->Add(fh1PtHard);
     fHistList->Add(fh1PtHard_NoW);
     fHistList->Add(fh1PtHard_Trials);
@@ -345,7 +404,10 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
     if (h1){
       // Printf("%s ",h1->GetName());
       h1->Sumw2();
+      continue;
     }
+    THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
+    if(hn)hn->Sumw2();
   }
 
   TH1::AddDirectory(oldStatus);
index 47057dae35da0aedca11d48ed3f243aeefdab964..bcf4e52ff8abf2328b6067d34760de0a11ac07d0 100644 (file)
@@ -5,8 +5,7 @@
  * See cxx source for full Copyright notice                               */
  
 #include "AliAnalysisTaskSE.h"
-////////////////
-#include <THnSparse.h>
+#include "THnSparse.h"
 ////////////////
 class AliJetHeader;
 class AliESDEvent;
@@ -18,6 +17,8 @@ class TList;
 class TChain;
 class TH2F;
 class TH3F;
+class TProfile;
+
 
 class AliAnalysisTaskJetSpectrum : public AliAnalysisTaskSE
 {
@@ -31,6 +32,8 @@ class AliAnalysisTaskJetSpectrum : public AliAnalysisTaskSE
     virtual void LocalInit() { Init(); }
     virtual void UserExec(Option_t *option);
     virtual void Terminate(Option_t *option);
+    virtual Bool_t Notify();
+
 
     virtual void SetExternalWeight(Float_t f){fExternalWeight = f;}
     virtual void SetUseExternalWeightOnly(Bool_t b){fUseExternalWeightOnly = b;}
@@ -74,6 +77,8 @@ class AliAnalysisTaskJetSpectrum : public AliAnalysisTaskSE
     Int_t         fAnalysisType;
     Float_t       fExternalWeight;    
 
+    TProfile*     fh1Xsec;   // pythia cross section and trials
+    TH1F*         fh1Trials; // trials are added
     TH1F*         fh1PtHard;  // Pt har of the event...       
     TH1F*         fh1PtHard_NoW;  // Pt har of the event...       
     TH1F*         fh1PtHard_Trials;  // Number of trials 
index 6f577464c95f498fca39bf42db49f524730fab8b..ebe4b66d1858536d194196b86f47f1d8e5e30144 100644 (file)
@@ -86,6 +86,8 @@ void unfold(const char* fileNameGen = "gen_pwg4spec.root", const char* folder =
   TFile::Open(fileNameRec);
   TH2F* hist = (TH2F*) gFile->Get("unfolding/fRecSpectrum");
   jetSpec->SetRecSpectrum(hist);
+  hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
+  if(hist->GetEntries()>0)jetSpec->SetGenSpectrum(hist);
 
   jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0);
   // last parameter = calculateErrors  <- this method to calculate the errors takes a lot of time
@@ -400,6 +402,24 @@ void FillSpecFromFiles(const char* fileNameReal = "histos_pwg4spec.root",const c
        jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
     }
 
+
+  // for control again, but now from the rec file
+  // generated jets (true distribution)
+  TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen"));  
+  for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++)
+    for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++)
+    {
+       Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te);
+       Float_t  z = fhEGenZGen->GetYaxis()->GetBinCenter(tz);
+       Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej);
+       Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z);
+       Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz);
+       Float_t err  = jetSpec->GetGenSpectrum()->GetBinError(bine,binz);
+       jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
+       jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
+    }
+
+
   file = TFile::Open("rec_pwg4spec.root", "RECREATE");
   jetSpec->SaveHistograms();
   file->Close();
@@ -409,13 +429,12 @@ void FillSpecFromFiles(const char* fileNameReal = "histos_pwg4spec.root",const c
   
 }
 
-
-
-
 void correct(){
   // simple steering to correct a given distribution;
   loadlibs();
-  FillSpecFromFiles("pwg4spec_0000000-0010000.root","pwg4spec_15-50.root");
+  // rec and gen
+  //  FillSpecFromFiles("pwg4spec_15-50_all.root","pwg4spec_allpt.root");
+  FillSpecFromFiles("pwg4spec_allpt.root","pwg4spec_allpt.root");
 
   char name[100];
   sprintf(name, "unfolded_pwg4spec.root");
@@ -424,3 +443,78 @@ void correct(){
   //draw(name, "unfolding", 1); 
 
 }
+
+void mergeJetAnaOutput(){
+  // This is used to merge the analysis-output from different 
+  // data samples/pt_hard bins
+  // in case the eventweigth was set to xsection/ntrials already this
+  // is not needed. Both methods only work in cse we do not mix different 
+  // pt_hard bins, and do not have overlapping bins
+
+  const Int_t nBins = 2;
+  // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05
+  // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03
+  // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02
+  const Float_t xsection[nBins] = {2.168291,2.44068e-02};
+  Float_t nTrials[nBins] = {0,0};
+  Float_t sf[nBins] = {0,0};
+
+  const char *cFile[nBins] = {"pwg4spec_15-50_all.root","pwg4spec_50_all.root"};
+
+
+  TList *lIn[nBins];
+  TFile *fIn[nBins];
+  for(int ib = 0;ib < nBins;++ib){
+    fIn[ib] = TFile::Open(cFile[ib]);
+    lIn[ib] = (TList*)fIn[ib]->Get("pwg4spec");
+    TH1* hTrials = (TH1F*)lIn[ib]->FindObject("fh1PtHard_Trials");
+    nTrials[ib] = hTrials->Integral();
+    sf[ib] = xsection[ib]/ nTrials[ib];
+  }
+
+  TFile *fOut = new TFile("pwg4spec_allpt.root","RECREATE");
+  TList *lOut = new TList();
+  lOut->SetName(lIn[0]->GetName());
+  // for the start scale all...
+  for(int ie = 0; ie < lIn[0]->GetEntries();++ie){
+    TH1 *h1Add = 0;
+    THnSparse *hnAdd = 0;
+    Printf("%d: %s",ie, lIn[0]->At(ie)->GetName());
+    for(int ib = 0;ib < nBins;++ib){
+      // dynamic cast does not work with cint
+      TObject *h = lIn[ib]->At(ie);
+      if(h->InheritsFrom("TH1")){
+       Printf("ib %d",ib);
+       TH1 *h1 = (TH1*)h;
+       if(ib==0){
+         h1Add = (TH1*)h1->Clone(h1->GetName());
+         h1Add->Scale(sf[ib]);
+       }
+       else{
+         h1Add->Add(h1,sf[ib]);
+       }
+      }
+      else if(h->InheritsFrom("THnSparse")){
+       Printf("ib %d",ib);
+       THnSparse *hn = (THnSparse*)h;
+       if(ib==0){
+         hnAdd = (THnSparse*)hn->Clone(hn->GetName());
+         hnAdd->Scale(sf[ib]);
+       }
+       else{
+         hnAdd->Add(hn,sf[ib]);
+       }
+      }
+      
+
+    }// ib
+    if(h1Add)lOut->Add(h1Add);
+    else if(hnAdd)lOut->Add(hnAdd);
+  }
+  fOut->cd();
+  lOut->Write(lOut->GetName(),TObject::kSingleKey);
+  fOut->Close();
+
+
+
+}