]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/JetTasks/AliAnalysisTaskJFSystematics.cxx
New analysis devoted to shower shape studies
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJFSystematics.cxx
index f07a9034600fd5376da526e9e4c04dea5afbd71c..90fbd3ca5044bb5a32f0356053185253483c6b7c 100644 (file)
@@ -79,6 +79,7 @@ AliAnalysisTaskJFSystematics::AliAnalysisTaskJFSystematics(): AliAnalysisTaskSE(
   fAnalysisType(0),
   fExternalWeight(1),    
   fRecEtaWindow(0.5),
+  fAvgTrials(1),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
@@ -126,6 +127,7 @@ AliAnalysisTaskJFSystematics::AliAnalysisTaskJFSystematics(const char* name):
   fAnalysisType(0),
   fExternalWeight(1),    
   fRecEtaWindow(0.5),
+  fAvgTrials(1),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
@@ -171,10 +173,12 @@ Bool_t AliAnalysisTaskJFSystematics::Notify()
   // Implemented Notify() to read the cross sections
   // and number of trials from pyxsec.root
   // 
+
+  fAvgTrials = 1; // reset for each file
+
   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
-  Double_t xsection = 0;
-  UInt_t   ntrials  = 0;
-  Float_t   ftrials  = 0;
+  Float_t xsection = 0;
+  Float_t ftrials  = 1;
   if(tree){
     TFile *curfile = tree->GetCurrentFile();
     if (!curfile) {
@@ -185,59 +189,11 @@ Bool_t AliAnalysisTaskJFSystematics::Notify()
       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
       return kFALSE;
     }
-
-    TString fileName(curfile->GetName());
-    if(fileName.Contains("AliESDs.root")){
-        fileName.ReplaceAll("AliESDs.root", "");
-    }
-    else if(fileName.Contains("AliAOD.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", "");
-    }
-    TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"));
-    if(!fxsec){
-      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);
-      }
-    }
-    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);
-    }
+    AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
     fh1Xsec->Fill("<#sigma>",xsection);
-    fh1Trials->Fill("#sum{ntrials}",ftrials);
+    // construct a poor man average trials 
+    Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
+    if(ftrials>=nEntries)fAvgTrials = ftrials/nEntries; 
   }
   return kTRUE;
 }
@@ -249,38 +205,6 @@ void AliAnalysisTaskJFSystematics::UserCreateOutputObjects()
   // Create the output container
   //
 
-  
-  // Connect the AOD
-
-  if(fUseAODInput){
-    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
-    if(!fAOD){
-      Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
-      return;
-    }
-    // fetch the header
-    fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
-    if(!fJetHeaderRec){
-      Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
-    }
-  }
-  else{
-    //  assume that the AOD is in the general output...
-    fAOD  = AODEvent();
-    if(!fAOD){
-      Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
-      return;
-    }
-    fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));    
-    if(!fJetHeaderRec){
-      Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
-    }
-    else{
-      if(fDebug>10)fJetHeaderRec->Dump();
-    }
-  }
-
 
   if (fDebug > 1) printf("AnalysisTaskJFSystematics::UserCreateOutputObjects() \n");
 
@@ -301,7 +225,7 @@ void AliAnalysisTaskJFSystematics::UserCreateOutputObjects()
       binLimitsPt[iPt] = 0.0;
     }
     else {// 1.0
-      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
+      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
     }
   }
   
@@ -329,7 +253,7 @@ void AliAnalysisTaskJFSystematics::UserCreateOutputObjects()
   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 = new TH1F("fh1Trials","trials event header or pyxsec file",1,0,1);
   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 
   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
@@ -386,19 +310,21 @@ void AliAnalysisTaskJFSystematics::UserCreateOutputObjects()
   fHistList->Add(fh1NRecJets);
   fHistList->Add(fh1PtRecIn);
   fHistList->Add(fh1PtRecOut);
-  fHistList->Add(fh1PtGenIn);
-  fHistList->Add(fh1PtGenOut);
-  fHistList->Add(fh2PtFGen);
-  fHistList->Add(fh2PhiFGen);
-  fHistList->Add(fh2EtaFGen);
-  fHistList->Add(fh2PtGenDeltaEta);
-  fHistList->Add(fh2PtGenDeltaPhi);
-  fHistList->Add(fh3RecOutEtaPhiPt);
-  fHistList->Add(fh3GenOutEtaPhiPt);      
-  fHistList->Add(fh3RecInEtaPhiPt);
-  fHistList->Add(fh3GenInEtaPhiPt);
-  fHistList->Add(fhnCorrelation);
 
+  if(fBranchGen.Length()>0){
+    fHistList->Add(fh1PtGenIn);
+    fHistList->Add(fh1PtGenOut);
+    fHistList->Add(fh2PtFGen);
+    fHistList->Add(fh2PhiFGen);
+    fHistList->Add(fh2EtaFGen);
+    fHistList->Add(fh2PtGenDeltaEta);
+    fHistList->Add(fh2PtGenDeltaPhi);
+    fHistList->Add(fh3RecOutEtaPhiPt);
+    fHistList->Add(fh3GenOutEtaPhiPt);      
+    fHistList->Add(fh3RecInEtaPhiPt);
+    fHistList->Add(fh3GenInEtaPhiPt);
+    fHistList->Add(fhnCorrelation);
+  }
 
   if(fAnalysisType==kSysJetOrder){
     // 
@@ -407,12 +333,15 @@ void AliAnalysisTaskJFSystematics::UserCreateOutputObjects()
       fHistList->Add(hTmp);
       hTmp = (TH1F*)fh1PtRecOut->Clone(Form("%s_%s%d",fh1PtRecOut->GetName(),fgkSysName[kSysJetOrder],i));
       fHistList->Add(hTmp);
-      hTmp = (TH1F*)fh1PtGenIn->Clone(Form("%s_%s%d",fh1PtGenIn->GetName(),fgkSysName[kSysJetOrder],i));
-      fHistList->Add(hTmp);
-      hTmp = (TH1F*)fh1PtGenOut->Clone(Form("%s_%s%d",fh1PtGenOut->GetName(),fgkSysName[kSysJetOrder],i));
-      fHistList->Add(hTmp);
-      THnSparseF *hnTmp = (THnSparseF*)fhnCorrelation->Clone(Form("%s_%s%d",fhnCorrelation->GetName(),fgkSysName[kSysJetOrder],i));
-      fHistList->Add(hnTmp);
+
+      if(fBranchGen.Length()>0){
+       hTmp = (TH1F*)fh1PtGenIn->Clone(Form("%s_%s%d",fh1PtGenIn->GetName(),fgkSysName[kSysJetOrder],i));
+       fHistList->Add(hTmp);
+       hTmp = (TH1F*)fh1PtGenOut->Clone(Form("%s_%s%d",fh1PtGenOut->GetName(),fgkSysName[kSysJetOrder],i));
+       fHistList->Add(hTmp);
+       THnSparseF *hnTmp = (THnSparseF*)fhnCorrelation->Clone(Form("%s_%s%d",fhnCorrelation->GetName(),fgkSysName[kSysJetOrder],i));
+       fHistList->Add(hnTmp);
+      }
     }
   }
 
@@ -448,9 +377,24 @@ void AliAnalysisTaskJFSystematics::UserExec(Option_t */*option*/)
   //
   // Execute analysis for current event
   //
-
-
-
+  if(fUseAODInput){    
+    fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+    if(!fAOD){
+      Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
+      return;
+    }
+    // fethc the header
+  }
+  else{
+    //  assume that the AOD is in the general output...
+    fAOD  = AODEvent();
+    if(!fAOD){
+      Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
+      return;
+    }
+  }
+  
   if (fDebug > 1)printf("AliAnalysisTaskJFSystematics::Analysing event # %5d\n", (Int_t) fEntry);
 
   // ========= These pointers need to be valid in any case ======= 
@@ -522,6 +466,9 @@ void AliAnalysisTaskJFSystematics::UserExec(Option_t */*option*/)
     }
   }// if we had the MCEvent
 
+  
+  fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
+  
   fh1PtHard->Fill(ptHard,eventW);
   fh1PtHardNoW->Fill(ptHard,1);
   fh1PtHardTrials->Fill(ptHard,nTrials);
@@ -558,8 +505,7 @@ void AliAnalysisTaskJFSystematics::UserExec(Option_t */*option*/)
   // Fetch the reconstructed jets...
   //
 
-  nRecJets = aodRecJets->GetEntries();
-  fh1NRecJets->Fill(nRecJets);
+
   nRecJets = TMath::Min(nRecJets,kMaxJets);
 
   for(int ir = 0;ir < nRecJets;++ir){