]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx
New book-keeping class
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum.cxx
index 96924c91c6ef554824f812d5c85aa16a94f7520b..b273deb738dac8236f41c6ff9c5f444a16a07866 100644 (file)
@@ -26,6 +26,7 @@
 #include <TInterpreter.h>
 #include <TChain.h>
 #include <TFile.h>
+#include <TKey.h>
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TH3F.h>
@@ -42,7 +43,6 @@
 #include "AliJetReader.h"
 #include "AliJetReaderHeader.h"
 #include "AliUA1JetHeaderV1.h"
-#include "AliJet.h"
 #include "AliESDEvent.h"
 #include "AliAODEvent.h"
 #include "AliAODHandler.h"
@@ -56,6 +56,7 @@
 #include "AliGenCocktailEventHeader.h"
 #include "AliInputEventHandler.h"
 
+
 #include "AliAnalysisHelperJetTasks.h"
 
 ClassImp(AliAnalysisTaskJetSpectrum)
@@ -67,9 +68,7 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
   fJetHeaderGen(0x0),
   fAOD(0x0),
   fBranchRec("jets"),
-  fConfigRec("ConfigJets.C"),
   fBranchGen(""),
-  fConfigGen(""),
   fUseAODInput(kFALSE),
   fUseExternalWeightOnly(kFALSE),
   fLimitGenJetEta(kFALSE),
@@ -83,7 +82,7 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
   fh1PtHardTrials(0x0),
   fh1NGenJets(0x0),
   fh1NRecJets(0x0),
-  fHistList(0x0) , 
+  fHistList(0x0) ,
   ////////////////
   fh1JetMultiplicity(0x0) ,     
   fh2ERecZRec(0x0) ,
@@ -110,9 +109,7 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
   fJetHeaderGen(0x0),
   fAOD(0x0),
   fBranchRec("jets"),
-  fConfigRec("ConfigJets.C"),
   fBranchGen(""),
-  fConfigGen(""),
   fUseAODInput(kFALSE),
   fUseExternalWeightOnly(kFALSE),
   fLimitGenJetEta(kFALSE),
@@ -161,6 +158,7 @@ Bool_t AliAnalysisTaskJetSpectrum::Notify()
   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
   Double_t xsection = 0;
   UInt_t   ntrials  = 0;
+  Float_t   ftrials  = 0;
   if(tree){
     TFile *curfile = tree->GetCurrentFile();
     if (!curfile) {
@@ -174,30 +172,56 @@ Bool_t AliAnalysisTaskJetSpectrum::Notify()
 
     TString fileName(curfile->GetName());
     if(fileName.Contains("AliESDs.root")){
-        fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
+        fileName.ReplaceAll("AliESDs.root", "");
     }
     else if(fileName.Contains("AliAOD.root")){
-        fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
+        fileName.ReplaceAll("AliAOD.root", "");
+    }
+    else if(fileName.Contains("AliAODs.root")){
+        fileName.ReplaceAll("AliAODs.root", "");
     }
     else if(fileName.Contains("galice.root")){
         // for running with galice and kinematics alone...                      
-        fileName.ReplaceAll("galice.root", "pyxsec.root");
+        fileName.ReplaceAll("galice.root", "");
     }
-    TFile *fxsec = TFile::Open(fileName.Data());
+    TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"));
     if(!fxsec){
-      Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
-      // no a severe condition
-      return kTRUE;
+      if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec.root"));
+      // next trial fetch the histgram file
+      fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
+      if(!fxsec){
+       // not a severe condition
+       if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec_hists.root"));        
+       return kTRUE;
+      }
+      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){
+         if(fDebug>0)Printf("%s:%d key not found in the file",(char*)__FILE__,__LINE__);       
+         return kTRUE;
+       }
+       TList *list = dynamic_cast<TList*>(key->ReadObj());
+       if(!list){
+         if(fDebug>0)Printf("%s:%d key is not a tlist",(char*)__FILE__,__LINE__);      
+         return kTRUE;
+       }
+       xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
+       ftrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
+      }
     }
-    TTree *xtree = (TTree*)fxsec->Get("Xsection");
-    if(!xtree){
-      Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+    else{
+      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);
+      ftrials = ntrials;
+      xtree->GetEntry(0);
     }
-    xtree->SetBranchAddress("xsection",&xsection);
-    xtree->SetBranchAddress("ntrials",&ntrials);
-    xtree->GetEntry(0);
     fh1Xsec->Fill("<#sigma>",xsection);
-    fh1Trials->Fill("#sum{ntrials}",ntrials);
+    fh1Trials->Fill("#sum{ntrials}",ftrials);
   }
   
   return kTRUE;
@@ -489,20 +513,8 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
     return;
   }
 
-
-  //  aodH->SelectEvent(kTRUE);
-
-  // ========= These pointers need to be valid in any case ======= 
-
-
-  /*
-  AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
-  if(!jhRec){
-    Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
-    return;
-  }
-  */
   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+
   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
   if(!aodRecJets){
     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
@@ -656,7 +668,7 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
   }
 
 
-  GetClosestJets(genJets,nGenJets,recJets,nRecJets,
+  AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
                 iGenIndex,iRecIndex,fDebug);
   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 
@@ -810,118 +822,3 @@ void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
 //
     if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
 }
-
-
-void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,const Int_t &nGenJets,
-                                               AliAODJet *recJets,const Int_t &nRecJets,
-                                               Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
-
-  //
-  // Relate the two input jet Arrays
-  //
-
-  //
-  // The association has to be unique
-  // So check in two directions
-  // find the closest rec to a gen
-  // 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
-  //    1   2   3
-  // ------------
-  // a| 3   2   0
-  // b| 0   1   0
-  // c| 0   0   3
-  // d| 0   0   1
-  // e| 0   0   1
-  // Topology
-  //   1     2
-  //     a         b        
-  //
-  //  d      c
-  //        3     e
-  // Only entries with "3" match from both sides
-  
-  const int kMode = 3;
-
-  Int_t iFlag[kMaxJets][kMaxJets];
-
-
-
-  for(int i = 0;i < kMaxJets;++i){
-    iRecIndex[i] = -1;
-    iGenIndex[i] = -1;
-    for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
-  }
-
-  if(nRecJets==0)return;
-  if(nGenJets==0)return;
-
-  const Float_t maxDist = 0.5;
-  // find the closest distance to the generated
-  for(int ig = 0;ig<nGenJets;++ig){
-    Float_t dist = maxDist;
-    if(iDebug>1)Printf("Gen (%d) p_T %3.3f eta %3.3f ph %3.3f ",ig,genJets[ig].Pt(),genJets[ig].Eta(),genJets[ig].Phi());
-    for(int ir = 0;ir<nRecJets;++ir){
-      Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
-      if(iDebug>1)Printf("Rec (%d) p_T %3.3f eta %3.3f ph %3.3f ",ir,recJets[ir].Pt(),recJets[ir].Eta(),recJets[ir].Phi());
-      if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
-      if(dR<dist){
-       iRecIndex[ig] = ir;
-       dist = dR;
-      }        
-    }
-    if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
-    // reset...
-    iRecIndex[ig] = -1;
-  }
-  // other way around
-  for(int ir = 0;ir<nRecJets;++ir){
-    Float_t dist = maxDist;
-    for(int ig = 0;ig<nGenJets;++ig){
-      Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
-      if(dR<dist){
-       iGenIndex[ir] = ig;
-       dist = dR;
-      }        
-    }
-    if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
-    // reset...
-    iGenIndex[ir] = -1;
-  }
-
-  // check for "true" correlations
-
-  if(iDebug>1)Printf(">>>>>> Matrix");
-
-  for(int ig = 0;ig<nGenJets;++ig){
-    for(int ir = 0;ir<nRecJets;++ir){
-      // Print
-      if(iDebug>1)printf("XFL %d ",iFlag[ig][ir]);
-
-      if(kMode==3){
-       // we have a uniqie correlation
-       if(iFlag[ig][ir]==3){
-         iGenIndex[ir] = ig;
-         iRecIndex[ig] = ir;
-       }
-      }
-      else{
-       // we just take the correlation from on side
-       if((iFlag[ig][ir]&2)==2){
-         iGenIndex[ir] = ig;
-       }
-       if((iFlag[ig][ir]&1)==1){
-         iRecIndex[ig] = ir;
-       }
-      }
-    }
-    if(iDebug>1)printf("\n");
-  }
-}
-
-
-