[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetCore.cxx

#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;  
}
Commit	Line	Data
75bf77e3	1	#include "TChain.h"
	2	#include "TTree.h"
	3	#include "TMath.h"
	4	#include "TH1F.h"
	5	#include "TH2F.h"
	6	#include "TH3F.h"
	7	#include "THnSparse.h"
	8	#include "TCanvas.h"
	9
	10	#include "AliLog.h"
	11
	12	#include "AliAnalysisTask.h"
	13	#include "AliAnalysisManager.h"
	14
	15	#include "AliVEvent.h"
	16	#include "AliESDEvent.h"
	17	#include "AliESDInputHandler.h"
	18	#include "AliCentrality.h"
	19	#include "AliAnalysisHelperJetTasks.h"
	20	#include "AliInputEventHandler.h"
	21	#include "AliAODJetEventBackground.h"
	22	#include "AliAnalysisTaskFastEmbedding.h"
	23
	24	#include "AliAODEvent.h"
	25	#include "AliAODJet.h"
	26
	27	#include "AliAnalysisTaskJetCore.h"
	28
	29	ClassImp(AliAnalysisTaskJetCore)
	30
	31	AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
	32	AliAnalysisTaskSE(),
	33	fESD(0x0),
	34	fAOD(0x0),
	35	fBackgroundBranch(""),
	36	fIsPbPb(kTRUE),
	37	fOfflineTrgMask(AliVEvent::kAny),
	38	fMinContribVtx(1),
	39	fVtxZMin(-8.),
	40	fVtxZMax(8.),
	41	fEvtClassMin(0),
	42	fEvtClassMax(4),
	43	fCentMin(0.),
	44	fCentMax(100.),
	45	fNInputTracksMin(0),
	46	fNInputTracksMax(-1),
	47	fJetEtaMin(-.5),
	48	fJetEtaMax(.5),
	49	fJetPtMin(20.),
	50	fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
	51	fJetPtFractionMin(0.5),
	52	fNMatchJets(4),
	53	fMatchMaxDist(0.8),
	54	fKeepJets(kFALSE),
	55	fRadioFrac(0.2),
	56	fMinDist(0.1),
	57	fkNbranches(2),
	58	fkEvtClasses(12),
	59	fOutputList(0x0),
	60	fbEvent(kTRUE),
	61	fHistEvtSelection(0x0),
	62	fHistJetSelection(0x0),
	63	fh2JetSelection(0x0),
	64	fh2JetCoreMethod1C10(0x0),
65	fh2JetCoreMethod2C10(0x0),
66	fh2JetCoreMethod3C10(0x0),
67	fh2JetCoreMethod1C30(0x0),
68	fh2JetCoreMethod2C30(0x0),
69	fh2JetCoreMethod3C30(0x0),
70	fh2JetCoreMethod1C60(0x0),
71	fh2JetCoreMethod2C60(0x0),
72	fh2JetCoreMethod3C60(0x0)
73
74
75	{
76	// default Constructor
77
78	fJetBranchName[0] = "";
79	fJetBranchName[1] = "";
80
81	fListJets[0] = new TList;
82	fListJets[1] = new TList;
83	}
84
85	AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
86	AliAnalysisTaskSE(name),
87	fESD(0x0),
88	fAOD(0x0),
89	fBackgroundBranch(""),
90	fIsPbPb(kTRUE),
91	fOfflineTrgMask(AliVEvent::kAny),
92	fMinContribVtx(1),
93	fVtxZMin(-8.),
94	fVtxZMax(8.),
95	fEvtClassMin(0),
96	fEvtClassMax(4),
97	fCentMin(0.),
98	fCentMax(100.),
99	fNInputTracksMin(0),
100	fNInputTracksMax(-1),
101	fJetEtaMin(-.5),
102	fJetEtaMax(.5),
103	fJetPtMin(20.),
104	fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
105	fJetPtFractionMin(0.5),
106	fNMatchJets(4),
107	fMatchMaxDist(0.8),
108	fKeepJets(kFALSE),
109	fRadioFrac(0.2),
110	fMinDist(0.1),
111	fkNbranches(2),
112	fkEvtClasses(12),
113	fOutputList(0x0),
114	fbEvent(kTRUE),
115	fHistEvtSelection(0x0),
116	fHistJetSelection(0x0),
117	fh2JetSelection(0x0),
118	fh2JetCoreMethod1C10(0x0),
119	fh2JetCoreMethod2C10(0x0),
120	fh2JetCoreMethod3C10(0x0),
121	fh2JetCoreMethod1C30(0x0),
122	fh2JetCoreMethod2C30(0x0),
123	fh2JetCoreMethod3C30(0x0),
124	fh2JetCoreMethod1C60(0x0),
125	fh2JetCoreMethod2C60(0x0),
126	fh2JetCoreMethod3C60(0x0)
127
128
129	{
130	// Constructor
131
132	fJetBranchName[0] = "";
133	fJetBranchName[1] = "";
134
135	fListJets[0] = new TList;
136	fListJets[1] = new TList;
137
138	DefineOutput(1, TList::Class());
139	}
140
141	AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
142	{
143	delete fListJets[0];
144	delete fListJets[1];
145	}
146
147	void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
148	{
149	fJetBranchName[0] = branch1;
150	fJetBranchName[1] = branch2;
151	}
152
153	void AliAnalysisTaskJetCore::Init()
154	{
155
156	// check for jet branches
157	if(!strlen(fJetBranchName[0].Data()) \|\| !strlen(fJetBranchName[1].Data())){
158	AliError("Jet branch name not set.");
159	}
160
161	}
162
163	void AliAnalysisTaskJetCore::UserCreateOutputObjects()
164	{
165	// Create histograms
166	// Called once
167	OpenFile(1);
168	if(!fOutputList) fOutputList = new TList;
169	fOutputList->SetOwner(kTRUE);
170
171	Bool_t oldStatus = TH1::AddDirectoryStatus();
172	TH1::AddDirectory(kFALSE);
173
174
175	fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
176	fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
177	fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
178	fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
179	fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
180	fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
181	fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
182
183	fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
184	fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
185	fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
186	fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
187	fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
188	fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
189	fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
190	fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
191	fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
192
193	fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
194	fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
195	fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
196	fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
197	fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
198	fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
199	fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
200	fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
201	fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
202
203
204	// UInt_t entries = 0; // bit coded, see GetDimParams() below
205	// UInt_t opt = 0; // bit coded, default (0) or high resolution (1)
206
207	fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
208	fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
209	fh2JetCoreMethod3C10 = new TH2F("JetCoreMethod3C10","",150, 0., 150.,100, 0., 1.5);
210	fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
211	fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
212	fh2JetCoreMethod3C30 = new TH2F("JetCoreMethod3C30","",150, 0., 150.,100, 0., 1.5);
213	fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
214	fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);
215	fh2JetCoreMethod3C60 = new TH2F("JetCoreMethod3C60","",150, 0., 150.,100, 0., 1.5);
216
217
218	fOutputList->Add(fHistEvtSelection);
219	fOutputList->Add(fHistJetSelection);
220	fOutputList->Add(fh2JetSelection);
221
222
223	fOutputList->Add(fh2JetCoreMethod1C10);
224	fOutputList->Add(fh2JetCoreMethod2C10);
225	fOutputList->Add(fh2JetCoreMethod3C10);
226	fOutputList->Add(fh2JetCoreMethod1C30);
227	fOutputList->Add(fh2JetCoreMethod2C30);
228	fOutputList->Add(fh2JetCoreMethod3C30);
229	fOutputList->Add(fh2JetCoreMethod1C60);
230	fOutputList->Add(fh2JetCoreMethod2C60);
231	fOutputList->Add(fh2JetCoreMethod3C60);
232
233
234
235	// =========== Switch on Sumw2 for all histos ===========
236	for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
237	TH1 h1 = dynamic_cast<TH1>(fOutputList->At(i));
238	if (h1){
239	h1->Sumw2();
240	continue;
241	}
242
243	}
244	TH1::AddDirectory(oldStatus);
245
246	PostData(1, fOutputList);
247	}
248
249	void AliAnalysisTaskJetCore::UserExec(Option_t *)
250	{
251
252
253	if(!strlen(fJetBranchName[0].Data()) \|\| !strlen(fJetBranchName[1].Data())){
254	AliError("Jet branch name not set.");
255	return;
256	}
257
258	fESD=dynamic_cast<AliESDEvent*>(InputEvent());
259	if (!fESD) {
260	AliError("ESD not available");
261	fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
262	} else {
263	fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
264	}
265	if (!fAOD) {
266	AliError("AOD not available");
267	return;
268	}
269
270	// -- event selection --
271	fHistEvtSelection->Fill(1); // number of events before event selection
272
273	// physics selection
274	AliInputEventHandler* inputHandler = (AliInputEventHandler*)
275	((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
276	if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
277	if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
278	fHistEvtSelection->Fill(2);
279	PostData(1, fOutputList);
280	return;
281	}
282
283	// vertex selection
284	AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
285	Int_t nTracksPrim = primVtx->GetNContributors();
286	if ((nTracksPrim < fMinContribVtx) \|\|
287	(primVtx->GetZ() < fVtxZMin) \|\|
288	(primVtx->GetZ() > fVtxZMax) ){
289	if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
290	fHistEvtSelection->Fill(3);
291	PostData(1, fOutputList);
292	return;
293	}
294
295	// event class selection (from jet helper task)
296	Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
297	if(fDebug) Printf("Event class %d", eventClass);
298	if (eventClass < fEvtClassMin \|\| eventClass > fEvtClassMax){
299	fHistEvtSelection->Fill(4);
300	PostData(1, fOutputList);
301	return;
302	}
303
304	// centrality selection
305	AliCentrality *cent = 0x0;
306	Float_t centValue = 0.;
307	if(fESD) cent = fESD->GetCentrality();
308	if(cent) centValue = cent->GetCentralityPercentile("V0M");
309	if(fDebug) printf("centrality: %f\n", centValue);
310	if (centValue < fCentMin \|\| centValue > fCentMax){
311	fHistEvtSelection->Fill(4);
312	PostData(1, fOutputList);
313	return;
314	}
315
316
317	// multiplicity due to input tracks
318	Int_t nInputTracks = GetNInputTracks();
319
320	if (nInputTracks < fNInputTracksMin \|\| (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
321	fHistEvtSelection->Fill(5);
322	PostData(1, fOutputList);
323	return;
324	}
325
326
327	fHistEvtSelection->Fill(0); // accepted events
328	// -- end event selection --
329
330
331	// get background
332	AliAODJetEventBackground* externalBackground = 0;
333	if(!externalBackground&&fBackgroundBranch.Length()){
334	externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
335	if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
336	}
337	Float_t rho = 0;
338	if(externalBackground)rho = externalBackground->GetBackground(0);
339
340
341	// fetch jets
342	TClonesArray *aodJets[2];
343	aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); //
344	aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); //
345
346
347
348	for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
349	fListJets[iJetType]->Clear();
350	if (!aodJets[iJetType]) continue;
351
352	if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
353
354	for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
355	AliAODJet jet = dynamic_cast<AliAODJet>((*aodJets[iJetType])[iJet]);
356	if (jet) fListJets[iJetType]->Add(jet);
357	}
358	fListJets[iJetType]->Sort();
359	}
360
361	Double_t etabig=0;
362	Double_t ptbig=0;
363	Double_t areabig=0;
364	Double_t phibig=0.;
365	Double_t etasmall=0;
366	Double_t ptsmall=0;
367	Double_t areasmall=0;
368	Double_t distr=0.;
369	Double_t phismall=0.;
370	for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
371	AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
372	etabig = jetbig->Eta();
373	phibig = jetbig->Phi();
374	ptbig = jetbig->Pt();
375	if(ptbig==0) continue;
376	areabig = jetbig->EffectiveAreaCharged();
377	if((etabig<fJetEtaMin)\|\|(etabig>fJetEtaMax)) continue;
378	Double_t dismin=100.;
379	Double_t ptmax=-10.;
380	Int_t index1=-1;
381	Int_t index2=-1;
382	Double_t fracin=0.;
383	//compute sum of the pt of the tracks in a concentric cone
384	TRefArray *genTrackList = jetbig->GetRefTracks();
385	Int_t nTracksGenJet = genTrackList->GetEntriesFast();
386	AliAODTrack* genTrack;
387	for(Int_t ir=0; ir<nTracksGenJet; ++ir){
388	genTrack = (AliAODTrack*)(genTrackList->At(ir));
389	Float_t etr=genTrack->Eta();
390	Float_t phir=genTrack->Phi();
391	distr=(etr-etabig)(etr-etabig)+(phir-phibig)(phir-phibig);
392	distr=TMath::Sqrt(distr);
393	if(distr<=fRadioFrac){ fracin=fracin+genTrack->Pt();}}
394	if(centValue<10) fh2JetCoreMethod3C10->Fill(ptbig-rho*areabig,fracin/ptbig);
395	if((centValue>30)&&(centValue<60)) fh2JetCoreMethod3C30->Fill(ptbig-rho*areabig,fracin/ptbig);
396	if(centValue>60) fh2JetCoreMethod3C60->Fill(ptbig-rho*areabig,fracin/ptbig);
397	//cout<<"method3 "<<fracin/ptbig<<" "<<ptbig-rho*areabig<<endl;
398	/////////////////////////////////////////////////////////////
399
400	for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
401	AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
402	etasmall = jetsmall->Eta();
403	phismall = jetsmall->Phi();
404	ptsmall = jetsmall->Pt();
405	areasmall = jetsmall->EffectiveAreaCharged();
406
407	if((etasmall<fJetEtaMin)\|\|(etasmall>fJetEtaMax)) continue;
408
409	Float_t tmpDeltaR=(phismall-phibig)(phismall-phibig)+(etasmall-etabig)(etasmall-etabig);
410	tmpDeltaR=TMath::Sqrt(tmpDeltaR);
411
412	if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
413	index2=j;}
414
415	if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
416	index1=j;}}
417
418	//method1:most concentric jet=core
419	if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
420
421	if(centValue<10) fh2JetCoreMethod1C10->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig);
422	if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig);
423	if(centValue>60) fh2JetCoreMethod1C60->Fill(ptbig-rho*areabig,jetmethod1->Pt()/ptbig);
424
425	//cout<<"method1 "<<jetmethod1->Pt()/ptbig<<" "<<ptbig<<endl;
426	}
427	//method2:hardest contained jet=core
428	if(index2!=-1){
429	AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
430	if(centValue<10) fh2JetCoreMethod2C10->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig);
431	if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig);
432	if(centValue>60) fh2JetCoreMethod2C60->Fill(ptbig-rho*areabig,jetmethod2->Pt()/ptbig);
433
434	//cout<<"method2 "<<jetmethod2->Pt()/ptbig<<" "<<ptbig<<endl;
435	}
436
437
438
439
440	}
441
442
443
444
445	PostData(1, fOutputList);
446	}
447
448	void AliAnalysisTaskJetCore::Terminate(const Option_t *)
449	{
450	// Draw result to the screen
451	// Called once at the end of the query
452
453	if (!GetOutputData(1))
454	return;
455	}
456
457	Int_t AliAnalysisTaskJetCore::GetNInputTracks()
458	{
459
460	Int_t nInputTracks = 0;
461
462	TString jbname(fJetBranchName[1]);
463	//needs complete event, use jets without background subtraction
464	for(Int_t i=1; i<=3; ++i){
465	if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
466	}
467	// use only HI event
468	if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
469	if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
470
471	if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
472	TClonesArray tmpAODjets = dynamic_cast<TClonesArray>(fAOD->FindListObject(jbname.Data()));
473	if(!tmpAODjets){
474	Printf("Jet branch %s not found", jbname.Data());
475	Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
476	return -1;
477	}
478
479	for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
480	AliAODJet jet = dynamic_cast<AliAODJet>((*tmpAODjets)[iJet]);
481	if(!jet) continue;
482	TRefArray *trackList = jet->GetRefTracks();
483	Int_t nTracks = trackList->GetEntriesFast();
484	nInputTracks += nTracks;
485	if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
486	}
487	if(fDebug) Printf("---> input tracks: %d", nInputTracks);
488
489	return nInputTracks;
490	}
491
492
493
494
495