New task for analysis of jet core (L. Cunqueiro)
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 2 Dec 2011 06:43:21 +0000 (06:43 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 2 Dec 2011 06:43:21 +0000 (06:43 +0000)
PWG4/CMakelibPWG4JetTasks.pkg
PWG4/JetTasks/AliAnalysisTaskJetCore.cxx [new file with mode: 0644]
PWG4/JetTasks/AliAnalysisTaskJetCore.h [new file with mode: 0644]
PWG4/PWG4JetTasksLinkDef.h
PWG4/macros/AddTaskJetCore.C [new file with mode: 0644]
PWG4/macros/AnalysisTrainPWG4Jets.C

index 5a42af9..010c587 100644 (file)
@@ -45,7 +45,8 @@ set ( SRCS
     JetTasks/AliAnalysisTaskJetChem.cxx 
     JetTasks/AliAnalysisTaskFragmentationFunction.cxx 
     JetTasks/AliAnalysisTaskMinijet.cxx 
-    JetTasks/AliAnalysisTaskQGSep.cxx 
+    JetTasks/AliAnalysisTaskQGSep.cxx     
+    JetTasks/AliAnalysisTaskJetCore.cxx        
     JetTasks/AliAnalysisTaskJetsTM.cxx 
     JetTasks/AliAnalysisTaskJetResponse.cxx 
     JetTasks/AliPWG4HighPtTrackQA.cxx 
diff --git a/PWG4/JetTasks/AliAnalysisTaskJetCore.cxx b/PWG4/JetTasks/AliAnalysisTaskJetCore.cxx
new file mode 100644 (file)
index 0000000..5ff619c
--- /dev/null
@@ -0,0 +1,495 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TMath.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "THnSparse.h"
+#include "TCanvas.h"
+
+#include "AliLog.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliVEvent.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliCentrality.h"
+#include "AliAnalysisHelperJetTasks.h"
+#include "AliInputEventHandler.h"
+#include "AliAODJetEventBackground.h"
+#include "AliAnalysisTaskFastEmbedding.h"
+
+#include "AliAODEvent.h"
+#include "AliAODJet.h"
+
+#include "AliAnalysisTaskJetCore.h"
+
+ClassImp(AliAnalysisTaskJetCore)
+
+AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
+AliAnalysisTaskSE(),
+fESD(0x0),
+fAOD(0x0),
+fBackgroundBranch(""),
+fIsPbPb(kTRUE),
+fOfflineTrgMask(AliVEvent::kAny),
+fMinContribVtx(1),
+fVtxZMin(-8.),
+fVtxZMax(8.),
+fEvtClassMin(0),
+fEvtClassMax(4),
+fCentMin(0.),
+fCentMax(100.),
+fNInputTracksMin(0),
+fNInputTracksMax(-1),
+fJetEtaMin(-.5),
+fJetEtaMax(.5),
+fJetPtMin(20.),
+fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
+fJetPtFractionMin(0.5),
+fNMatchJets(4),
+fMatchMaxDist(0.8),
+fKeepJets(kFALSE),
+fRadioFrac(0.2),
+fMinDist(0.1),
+fkNbranches(2),
+fkEvtClasses(12),
+fOutputList(0x0),
+fbEvent(kTRUE),
+fHistEvtSelection(0x0),
+fHistJetSelection(0x0),
+fh2JetSelection(0x0),
+fh2JetCoreMethod1C10(0x0),
+fh2JetCoreMethod2C10(0x0),
+fh2JetCoreMethod3C10(0x0),
+fh2JetCoreMethod1C30(0x0),
+fh2JetCoreMethod2C30(0x0),
+fh2JetCoreMethod3C30(0x0),
+fh2JetCoreMethod1C60(0x0),
+fh2JetCoreMethod2C60(0x0),
+fh2JetCoreMethod3C60(0x0)
+
+{
+   // default Constructor
+
+   fJetBranchName[0] = "";
+   fJetBranchName[1] = "";
+
+   fListJets[0] = new TList;
+   fListJets[1] = new TList;
+}
+
+AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
+AliAnalysisTaskSE(name),
+fESD(0x0),
+fAOD(0x0),
+fBackgroundBranch(""),
+fIsPbPb(kTRUE),
+fOfflineTrgMask(AliVEvent::kAny),
+fMinContribVtx(1),
+fVtxZMin(-8.),
+fVtxZMax(8.),
+fEvtClassMin(0),
+fEvtClassMax(4),
+fCentMin(0.),
+fCentMax(100.),
+fNInputTracksMin(0),
+fNInputTracksMax(-1),
+fJetEtaMin(-.5),
+fJetEtaMax(.5),
+fJetPtMin(20.),
+fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
+fJetPtFractionMin(0.5),
+fNMatchJets(4),
+fMatchMaxDist(0.8),
+fKeepJets(kFALSE),
+fRadioFrac(0.2),
+fMinDist(0.1),
+fkNbranches(2),
+fkEvtClasses(12),
+fOutputList(0x0),
+fbEvent(kTRUE),
+fHistEvtSelection(0x0),
+fHistJetSelection(0x0),
+fh2JetSelection(0x0),
+fh2JetCoreMethod1C10(0x0),
+fh2JetCoreMethod2C10(0x0),
+fh2JetCoreMethod3C10(0x0),
+fh2JetCoreMethod1C30(0x0),
+fh2JetCoreMethod2C30(0x0),
+fh2JetCoreMethod3C30(0x0),
+fh2JetCoreMethod1C60(0x0),
+fh2JetCoreMethod2C60(0x0),
+fh2JetCoreMethod3C60(0x0)
+
+
+ {
+   // Constructor
+
+   fJetBranchName[0] = "";
+   fJetBranchName[1] = "";
+
+   fListJets[0] = new TList;
+   fListJets[1] = new TList;
+
+   DefineOutput(1, TList::Class());
+}
+
+AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
+{
+   delete fListJets[0];
+   delete fListJets[1];
+}
+
+void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
+{
+   fJetBranchName[0] = branch1;
+   fJetBranchName[1] = branch2;
+}
+
+void AliAnalysisTaskJetCore::Init()
+{
+
+   // check for jet branches
+   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
+      AliError("Jet branch name not set.");
+   }
+
+}
+
+void AliAnalysisTaskJetCore::UserCreateOutputObjects()
+{
+   // Create histograms
+   // Called once
+   OpenFile(1);
+   if(!fOutputList) fOutputList = new TList;
+   fOutputList->SetOwner(kTRUE);
+
+   Bool_t oldStatus = TH1::AddDirectoryStatus();
+   TH1::AddDirectory(kFALSE);
+
+
+   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
+   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
+   fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
+
+   fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
+   fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
+   fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
+   fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
+   fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
+   fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
+   fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
+   fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
+   fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
+
+   fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
+   fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
+   fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
+   fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
+   fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
+   fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
+   fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
+   fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
+   fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
+
+
+   //   UInt_t entries = 0; // bit coded, see GetDimParams() below
+   //   UInt_t opt = 0;  // bit coded, default (0) or high resolution (1)
+
+    fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
+    fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
+    fh2JetCoreMethod3C10 = new TH2F("JetCoreMethod3C10","",150, 0., 150.,100, 0., 1.5); 
+    fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
+    fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
+    fh2JetCoreMethod3C30 = new TH2F("JetCoreMethod3C30","",150, 0., 150.,100, 0., 1.5); 
+    fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
+    fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);
+    fh2JetCoreMethod3C60 = new TH2F("JetCoreMethod3C60","",150, 0., 150.,100, 0., 1.5); 
+
+
+   fOutputList->Add(fHistEvtSelection);
+   fOutputList->Add(fHistJetSelection);
+   fOutputList->Add(fh2JetSelection);
+  
+  
+      fOutputList->Add(fh2JetCoreMethod1C10);
+      fOutputList->Add(fh2JetCoreMethod2C10);
+      fOutputList->Add(fh2JetCoreMethod3C10);
+      fOutputList->Add(fh2JetCoreMethod1C30);
+      fOutputList->Add(fh2JetCoreMethod2C30);
+      fOutputList->Add(fh2JetCoreMethod3C30);
+      fOutputList->Add(fh2JetCoreMethod1C60);
+      fOutputList->Add(fh2JetCoreMethod2C60);
+      fOutputList->Add(fh2JetCoreMethod3C60);
+     
+     
+
+   // =========== Switch on Sumw2 for all histos ===========
+   for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
+      TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
+      if (h1){
+         h1->Sumw2();
+         continue;
+      }
+    
+   }
+   TH1::AddDirectory(oldStatus);
+
+   PostData(1, fOutputList);
+}
+
+void AliAnalysisTaskJetCore::UserExec(Option_t *)
+{
+   
+
+   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
+      AliError("Jet branch name not set.");
+      return;
+   }
+
+   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
+   if (!fESD) {
+      AliError("ESD not available");
+      fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+   } else {
+      fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
+   }
+   if (!fAOD) {
+      AliError("AOD not available");
+      return;
+   }
+
+   // -- event selection --
+   fHistEvtSelection->Fill(1); // number of events before event selection
+
+   // physics selection
+   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
+   ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
+      if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
+      fHistEvtSelection->Fill(2);
+      PostData(1, fOutputList);
+      return;
+   }
+
+   // vertex selection
+   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
+   Int_t nTracksPrim = primVtx->GetNContributors();
+   if ((nTracksPrim < fMinContribVtx) ||
+         (primVtx->GetZ() < fVtxZMin) ||
+         (primVtx->GetZ() > fVtxZMax) ){
+      if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
+      fHistEvtSelection->Fill(3);
+      PostData(1, fOutputList);
+      return;
+   }
+
+   // event class selection (from jet helper task)
+   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
+   if(fDebug) Printf("Event class %d", eventClass);
+   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
+      fHistEvtSelection->Fill(4);
+      PostData(1, fOutputList);
+      return;
+   }
+
+   // centrality selection
+   AliCentrality *cent = 0x0;
+   Float_t centValue = 0.; 
+   if(fESD) cent = fESD->GetCentrality();
+   if(cent) centValue = cent->GetCentralityPercentile("V0M");
+   if(fDebug) printf("centrality: %f\n", centValue);
+   if (centValue < fCentMin || centValue > fCentMax){
+      fHistEvtSelection->Fill(4);
+      PostData(1, fOutputList);
+      return;
+   }
+
+
+   // multiplicity due to input tracks
+   Int_t nInputTracks = GetNInputTracks();
+
+   if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
+      fHistEvtSelection->Fill(5);
+      PostData(1, fOutputList);
+      return;
+   }
+
+   
+   fHistEvtSelection->Fill(0); // accepted events  
+   // -- end event selection --
+
+
+   // get background
+   AliAODJetEventBackground* externalBackground = 0;
+   if(!externalBackground&&fBackgroundBranch.Length()){
+      externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
+      if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
+   }
+   Float_t rho = 0;
+   if(externalBackground)rho = externalBackground->GetBackground(0);
+
+
+   // fetch jets
+   TClonesArray *aodJets[2];
+   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // 
+   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // 
+
+
+
+   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
+      fListJets[iJetType]->Clear();
+      if (!aodJets[iJetType]) continue;
+
+      if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
+      
+      for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
+         AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
+         if (jet) fListJets[iJetType]->Add(jet);
+      }
+      fListJets[iJetType]->Sort();
+   }
+   
+   Double_t etabig=0;
+   Double_t ptbig=0;
+   Double_t areabig=0;
+   Double_t phibig=0.;
+   Double_t etasmall=0;
+   Double_t ptsmall=0;
+   Double_t areasmall=0;
+   Double_t distr=0.;
+   Double_t phismall=0.;
+   for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
+           AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
+           etabig  = jetbig->Eta();
+           phibig  = jetbig->Phi();
+           ptbig   = jetbig->Pt();
+           if(ptbig==0) continue; 
+           areabig = jetbig->EffectiveAreaCharged();
+          if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
+                   Double_t dismin=100.;
+                   Double_t ptmax=-10.; 
+                   Int_t index1=-1;
+                   Int_t index2=-1;
+                   Double_t fracin=0.; 
+          //compute sum of the pt of the tracks in a concentric cone
+           TRefArray *genTrackList = jetbig->GetRefTracks();
+           Int_t nTracksGenJet = genTrackList->GetEntriesFast();
+           AliAODTrack* genTrack;
+             for(Int_t ir=0; ir<nTracksGenJet; ++ir){
+             genTrack = (AliAODTrack*)(genTrackList->At(ir));
+            Float_t etr=genTrack->Eta();
+             Float_t phir=genTrack->Phi();
+             distr=(etr-etabig)*(etr-etabig)+(phir-phibig)*(phir-phibig);
+             distr=TMath::Sqrt(distr);
+            if(distr<=fRadioFrac){ fracin=fracin+genTrack->Pt();}}    
+            if(centValue<10) fh2JetCoreMethod3C10->Fill(ptbig-rho*areabig,fracin/ptbig);
+             if((centValue>30)&&(centValue<60)) fh2JetCoreMethod3C30->Fill(ptbig-rho*areabig,fracin/ptbig);
+             if(centValue>60) fh2JetCoreMethod3C60->Fill(ptbig-rho*areabig,fracin/ptbig); 
+             //cout<<"method3  "<<fracin/ptbig<<" "<<ptbig-rho*areabig<<endl;
+            /////////////////////////////////////////////////////////////
+
+                for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
+                  AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
+                  etasmall  = jetsmall->Eta();
+                  phismall = jetsmall->Phi();
+                  ptsmall   = jetsmall->Pt();
+                  areasmall = jetsmall->EffectiveAreaCharged();
+
+                  if((etasmall<fJetEtaMin)||(etasmall>fJetEtaMax)) continue;
+                  
+                   Float_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
+                  tmpDeltaR=TMath::Sqrt(tmpDeltaR);
+
+                 if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
+                   index2=j;}  
+  
+                  if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
+                   index1=j;}}
+
+                //method1:most concentric jet=core 
+               if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));       
+       
+                 if(centValue<10) fh2JetCoreMethod1C10->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig);
+                 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig);
+                 if(centValue>60) fh2JetCoreMethod1C60->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig); 
+
+                  //cout<<"method1  "<<jetmethod1->Pt()/ptbig<<" "<<ptbig<<endl;   
+               }
+                //method2:hardest contained jet=core   
+               if(index2!=-1){ 
+                  AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
+                  if(centValue<10) fh2JetCoreMethod2C10->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig);
+                 if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig);
+                 if(centValue>60) fh2JetCoreMethod2C60->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig); 
+
+                 //cout<<"method2  "<<jetmethod2->Pt()/ptbig<<" "<<ptbig<<endl;
+                               }  
+
+   
+      
+       
+   }
+     
+      
+
+
+   PostData(1, fOutputList);
+}
+
+void AliAnalysisTaskJetCore::Terminate(const Option_t *)
+{
+   // Draw result to the screen
+   // Called once at the end of the query
+
+   if (!GetOutputData(1))
+   return;
+}
+
+Int_t AliAnalysisTaskJetCore::GetNInputTracks()
+{
+
+   Int_t nInputTracks = 0;
+
+   TString jbname(fJetBranchName[1]);
+   //needs complete event, use jets without background subtraction
+   for(Int_t i=1; i<=3; ++i){
+      if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
+   }
+   // use only HI event
+   if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
+   if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
+
+   if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
+   TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
+   if(!tmpAODjets){
+      Printf("Jet branch %s not found", jbname.Data());
+      Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
+      return -1;
+   }
+   
+   for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
+      AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
+      if(!jet) continue;
+      TRefArray *trackList = jet->GetRefTracks();
+      Int_t nTracks = trackList->GetEntriesFast();
+      nInputTracks += nTracks;
+      if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
+   }
+   if(fDebug) Printf("---> input tracks: %d", nInputTracks);
+
+   return nInputTracks;  
+}
+
+
+
+
+
diff --git a/PWG4/JetTasks/AliAnalysisTaskJetCore.h b/PWG4/JetTasks/AliAnalysisTaskJetCore.h
new file mode 100644 (file)
index 0000000..00b6bca
--- /dev/null
@@ -0,0 +1,139 @@
+#ifndef ALIANALYSISTASKJETCORE_H
+#define ALIANALYSISTASKJETCORE_H
+
+class TH1F;
+class TH2F;
+class TH3F;
+class THnSparse;
+class AliESDEvent;
+class AliAODEvent;
+
+#include "AliAnalysisTaskSE.h"
+#include "AliVEvent.h"
+
+class AliAnalysisTaskJetCore : public AliAnalysisTaskSE {
+public:
+   AliAnalysisTaskJetCore();
+   AliAnalysisTaskJetCore(const char *name);
+   virtual ~AliAnalysisTaskJetCore();
+
+   virtual void     LocalInit() {Init();}
+   virtual void     Init();
+   virtual void     UserCreateOutputObjects();
+   virtual void     UserExec(Option_t *option);
+   virtual void     Terminate(const Option_t*);
+
+   virtual Int_t      GetNInputTracks();
+
+  
+
+   virtual AliVEvent::EOfflineTriggerTypes GetOfflineTrgMask() const { return fOfflineTrgMask; }
+   virtual void     GetBranchNames(TString &branch1, TString &branch2) const { branch1 = fJetBranchName[0]; branch2 = fJetBranchName[1]; }
+   virtual Bool_t   GetIsPbPb() const { return fIsPbPb; }
+   virtual Int_t    GetMinContribVtx() const { return fMinContribVtx; };
+   virtual Float_t  GetVtxZMin() const { return fVtxZMin; }
+   virtual Float_t  GetVtxZMax() const { return fVtxZMax; }
+   virtual Int_t    GetEvtClassMin() const { return fEvtClassMin; }
+   virtual Int_t    GetEvtClassMax() const { return fEvtClassMax; }
+   virtual Float_t  GetCentMin() const { return fCentMin; }
+   virtual Float_t  GetCentMax() const { return fCentMax; }
+   virtual Int_t    GetNInputTracksMin() const { return fNInputTracksMin; }
+   virtual Int_t    GetNInputTracksMax() const { return fNInputTracksMax; } 
+   virtual Float_t  GetJetEtaMin() const { return fJetEtaMin; }
+   virtual Float_t  GetJetEtaMax() const { return fJetEtaMax; }
+   virtual Float_t  GetJetPtMin() const { return fJetPtMin; }
+   virtual Float_t  GetJetPtFractionMin() const { return fJetPtFractionMin; }
+   virtual Int_t    GetNMatchJets() const { return fNMatchJets; }
+   virtual void     SetBranchNames(const TString &branch1, const TString &branch2);
+   virtual void     SetBackgroundBranch(TString &branch) { fBackgroundBranch = branch;}
+   virtual void     SetIsPbPb(Bool_t b=kTRUE) { fIsPbPb = b; }
+   virtual void     SetOfflineTrgMask(AliVEvent::EOfflineTriggerTypes mask) { fOfflineTrgMask = mask; }
+   virtual void     SetMinContribVtx(Int_t n) { fMinContribVtx = n; }
+   virtual void     SetVtxZMin(Float_t z) { fVtxZMin = z; }
+   virtual void     SetVtxZMax(Float_t z) { fVtxZMax = z; }
+   virtual void     SetEvtClassMin(Int_t evtClass) { fEvtClassMin = evtClass; }
+   virtual void     SetEvtClassMax(Int_t evtClass) { fEvtClassMax = evtClass; }
+   virtual void     SetCentMin(Float_t cent) { fCentMin = cent; }
+   virtual void     SetCentMax(Float_t cent) { fCentMax = cent; }
+   virtual void     SetNInputTracksMin(Int_t nTr) { fNInputTracksMin = nTr; }
+   virtual void     SetNInputTracksMax(Int_t nTr) { fNInputTracksMax = nTr; }
+   virtual void     SetJetEtaMin(Float_t eta) { fJetEtaMin = eta; }
+   virtual void     SetJetEtaMax(Float_t eta) { fJetEtaMax = eta; }
+   virtual void     SetJetPtMin(Float_t pt) { fJetPtMin = pt; }
+   virtual void     SetJetTriggerExclude(UChar_t i) { fJetTriggerExcludeMask = i; }
+   virtual void     SetJetPtFractionMin(Float_t frac) { fJetPtFractionMin = frac; }
+   virtual void     SetNMatchJets(Int_t n) { fNMatchJets = n; }
+   virtual void     SetFillEvent(Bool_t b) { fbEvent = b; }
+   virtual void     SetKeepJets(Bool_t b = kTRUE) { fKeepJets = b; }
+   virtual void     SetRadioFrac(Float_t radiofrac) { fRadioFrac = radiofrac; }
+   virtual void     SetMinDist(Float_t minDist) { fMinDist = minDist; }
+
+
+private:
+   // ESD/AOD events
+   AliESDEvent *fESD;    //! ESD object
+   AliAODEvent *fAOD;    //! AOD event
+
+   // jets to compare
+   TString fJetBranchName[2]; //  name of jet branches to compare
+   TList *fListJets[2];       //! jet lists
+
+   TString fBackgroundBranch;
+
+   // event selection
+   Bool_t fIsPbPb;         // is Pb-Pb (fast embedding) or p-p (detector response)
+   AliVEvent::EOfflineTriggerTypes fOfflineTrgMask; // mask of offline triggers to accept
+   Int_t   fMinContribVtx; // minimum number of track contributors for primary vertex
+   Float_t fVtxZMin;     // lower bound on vertex z
+   Float_t fVtxZMax;     // upper bound on vertex z
+   Int_t   fEvtClassMin;         // lower bound on event class
+   Int_t   fEvtClassMax;         // upper bound on event class
+   Float_t fCentMin;     // lower bound on centrality
+   Float_t fCentMax;     // upper bound on centrality
+   Int_t   fNInputTracksMin;  // lower bound of nb. of input tracks
+   Int_t   fNInputTracksMax;  // upper bound of nb. of input tracks
+   Float_t fJetEtaMin;     // lower bound on eta for found jets
+   Float_t fJetEtaMax;     // upper bound on eta for found jets
+   Float_t fJetPtMin;      // minimum jet pT
+   UChar_t fJetTriggerExcludeMask; // mask for jet triggeres to exclude
+   Float_t fJetPtFractionMin; // minimum fraction for positiv match of jets
+   Int_t   fNMatchJets;       // maximal nb. of jets taken for matching
+   Double_t fMatchMaxDist;     // maximal distance of matching jets
+   Bool_t  fKeepJets;          // keep jets with negative pt after background subtraction
+   Float_t fRadioFrac;                          //size of the concentric cone
+   Float_t fMinDist;                            //minimum jet axis distance
+
+   // output objects
+   const Int_t fkNbranches;                   //! number of branches to be read
+   const Int_t fkEvtClasses;                  //! number of event classes
+   
+   TList *fOutputList;                        //! output data container
+   Bool_t fbEvent;                            // fill fhnEvent
+   TH1I  *fHistEvtSelection;                  //! event selection statistic
+   TH1I  *fHistJetSelection;                  //! jet selection statistic
+   TH2F  *fh2JetSelection;                    //! jet selection statistic, with 
+  
+   TH2F      *fh2JetCoreMethod1C10;           //jet core method1 0-10%
+   TH2F      *fh2JetCoreMethod2C10;           //jet core method2 0-10%
+   TH2F      *fh2JetCoreMethod3C10;           //jet core method3 0-10%
+
+   TH2F      *fh2JetCoreMethod1C30;           //jet core method1 30-60%
+   TH2F      *fh2JetCoreMethod2C30;           //jet core method2 30-60%
+   TH2F      *fh2JetCoreMethod3C30;           //jet core method3 30-60%
+      
+
+   TH2F      *fh2JetCoreMethod1C60;          //jet core method1 60-100% 
+   TH2F      *fh2JetCoreMethod2C60;          //jet core method2 60-100%
+   TH2F      *fh2JetCoreMethod3C60;          //jet core method3 60-100%
+   
+
+
+
+   AliAnalysisTaskJetCore(const AliAnalysisTaskJetCore&); // not implemented
+   AliAnalysisTaskJetCore& operator=(const AliAnalysisTaskJetCore&); // not implemented
+
+   ClassDef(AliAnalysisTaskJetCore, 4);
+};
+
+#endif
+
index bc5381f..b5db74e 100644 (file)
@@ -32,6 +32,7 @@
 #pragma link C++ class AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass+;
 #pragma link C++ class AliAnalysisTaskMinijet+;
 #pragma link C++ class AliAnalysisTaskJetsTM+;
+#pragma link C++ class AliAnalysisTaskJetCore+;
 #pragma link C++ class AliAnalysisTaskQGSep+;
 #pragma link C++ class AliAnalysisTaskJetResponse+;
 #pragma link C++ class AliAnalysisTaskJetResponseV2;
diff --git a/PWG4/macros/AddTaskJetCore.C b/PWG4/macros/AddTaskJetCore.C
new file mode 100644 (file)
index 0000000..01ec490
--- /dev/null
@@ -0,0 +1,65 @@
+
+
+AliAnalysisTaskJetCore* AddTaskJetCore(const char* bRec1,const char* bRec2, UInt_t filterMask = 272 , Float_t ptTrackMin = 0.15,  Int_t eventClassMin = 0, Int_t eventClassMax = 4){
+
+   Printf("adding task jet core\n");
+
+   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+   if(!mgr){
+      ::Error("AddTaskJetCore", "No analysis manager to connect to.");
+      return NULL;
+   }
+   if(!mgr->GetInputEventHandler()){
+      ::Error("AddTaskJetCore", "This task requires an input event handler.");
+      return NULL;
+   }
+
+     
+   //jetsin bRec1 branch have larger radius than those in bRec2  
+
+  TString typeRec(bRec1);
+  TString typeGen(bRec2);
+      
+  AliAnalysisTaskJetCore *task = new AliAnalysisTaskJetCore(Form("JetCore_%s_%s",bRec1,bRec2));
+   
+
+
+   task->SetBranchNames(bRec1,bRec2);
+   task->SetOfflineTrgMask(AliVEvent::kMB);
+
+   task->SetEvtClassMin(eventClassMin);
+   task->SetEvtClassMax(eventClassMax);
+   task->SetCentMin(0.);
+   task->SetCentMax(100.);
+
+   
+   task->SetJetPtMin(0.);   
+   
+   //task->SetRadioFrac(0.2);   //radius of the concentric cone
+                                //within which you sum up the pT of
+                                //the tracks to compute the core fraction
+                                //in method3. It is also the maxium distance
+                                //between jet axis in method2. Default is 0.2 
+
+   //task->SetMinDist(0.1);     //minimum distance between jets to be 
+                                //concentric in method1. Default is 0.1
+
+   task->SetBackgroundBranch("jeteventbackground_clustersAOD_KT04_B0_Filter00272_Cut00150_Skip00"); 
+
+   mgr->AddTask(task);
+
+
+   AliAnalysisDataContainer *coutputJetCore = mgr->CreateContainer(Form("pwg4jetcore_%s_%s",bRec1,bRec2), TList::Class(),AliAnalysisManager::kOutputContainer,Form("%s:PWG4_jetcore_%s_%s",AliAnalysisManager::GetCommonFileName(),bRec1,bRec2));
+
+
+
+
+
+   mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer());
+   mgr->ConnectOutput(task, 0, mgr->GetCommonOutputContainer());
+   mgr->ConnectOutput(task, 1, coutputJetCore);
+
+   return task;
+}
index 65ef02b..9671f26 100644 (file)
@@ -911,7 +911,7 @@ void AnalysisTrainPWG4Jets(const char *analysis_mode="local",
      taskjetServ->SetPhysicsSelectionFlag(iPhysicsSelectionFlag); // 
      taskjetServ->SetNonStdFile(kDeltaAODJetName.Data());
      taskjetServ->SetTrackEtaWindow(fTrackEtaWindow);
-       taskjetServ->SetZVertexCut(fVertexWindow);
+     taskjetServ->SetZVertexCut(fVertexWindow);
      taskjetServ->SetFilterMask(kHighPtFilterMask);
 
      if(kIsPbPb)taskjetServ->SetCollisionType(AliAnalysisTaskJetServices::kPbPb);