From: kleinb Date: Tue, 22 Jul 2008 11:11:54 +0000 (+0000) Subject: Initial structure for PWG4JetTasks (for PWG4 Task Force) incl. steering macro and... X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=e989c76a0ca1f598e3b970fd2a1541f2c6f75e57 Initial structure for PWG4JetTasks (for PWG4 Task Force) incl. steering macro and modified Makefile. First Task: Underlying Event studies --- diff --git a/PWG4/AliAnalysisTaskUE.cxx b/PWG4/AliAnalysisTaskUE.cxx new file mode 100644 index 00000000000..7ca9cc504e7 --- /dev/null +++ b/PWG4/AliAnalysisTaskUE.cxx @@ -0,0 +1,725 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id:$ */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliAnalysisTaskUE.h" +#include "AliAnalysisManager.h" +#include "AliMCEventHandler.h" +#include "AliAODEvent.h" +#include "AliAODInputHandler.h" +#include "AliAODHandler.h" +#include "AliStack.h" +#include "AliAODJet.h" +#include "AliAODTrack.h" + +#include "AliLog.h" + +// +// Analysis class for Underlying Event studies +// +// Look for correlations on the tranverse regions to +// the leading charged jet +// +// This class needs as input AOD with track and Jets +// the output is a list of histograms +// +// AOD can be either connected to the InputEventHandler +// for a chain of AOD files +// or +// to the OutputEventHandler +// for a chain of ESD files, so this case class should be +// in the train after the Jet finder +// +// Arian.Abrahantes.Quintana@cern.ch +// Ernesto.Lopez.Torres@cern.ch +// + + +ClassImp( AliAnalysisTaskUE) + +//////////////////////////////////////////////////////////////////////// + + +//____________________________________________________________________ +AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name): + AliAnalysisTask(name, ""), + fDebug(kFALSE), + fAOD(0x0), + fListOfHistos(0x0), + fBinsPtInHist(30), + fMinJetPtInHist(0.), + fMaxJetPtInHist(300.), + fAnaType(1), + fRegionType(1), + fConeRadius(0.7), + fOrdering(1), + fJet1EtaCut(0.2), + fJet2DeltaPhiCut(2.616), // 150 degrees + fJet2RatioPtCut(0.8), + fJet3PtCut(15.), + fTrackPtCut(0.), + fTrackEtaCut(0.9), + fhNJets(0x0), + fhEleadingPt(0x0), + fhMinRegPtDist(0x0), + fhRegionMultMin(0x0), + fhMinRegAvePt(0x0), + fhMinRegSumPt(0x0), + fhMinRegMaxPtPart(0x0), + fhMinRegSumPtvsMult(0x0), + fhdNdEta_PhiDist(0x0), + fhFullRegPartPtDistVsEt(0x0), + fhTransRegPartPtDistVsEt(0x0), + fhRegionSumPtMaxVsEt(0x0), + fhRegionMultMax(0x0), + fhRegionMultMaxVsEt(0x0), + fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0), + fhRegionMultMinVsEt(0x0), + fhRegionAveSumPtVsEt(0x0), + fhRegionDiffSumPtVsEt(0x0), + fhRegionAvePartPtMaxVsEt(0x0), + fhRegionAvePartPtMinVsEt(0x0), + fhRegionMaxPartPtMaxVsEt(0x0)//, fhValidRegion(0x0) +{ + // Default constructor + // Define input and output slots here + // Input slot #0 works with a TChain + DefineInput(0, TChain::Class()); + // Output slot #0 writes into a TList container + DefineOutput(0, TList::Class()); +} + + +//____________________________________________________________________ +void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/) +{ +// Connect the input data + +// We need AOD with tracks and jets. +// Since AOD can be either connected to the InputEventHandler or to the OutputEventHandler +// we need to check where it is and get the pointer to AODEvent in the right way + + if (fDebug > 1) AliInfo("ConnectInputData() \n"); + + + TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); + + if( handler && handler->InheritsFrom("AliAODInputHandler") ) { + fAOD = ((AliAODInputHandler*)handler)->GetEvent(); + } else { + handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler(); + if( handler && handler->InheritsFrom("AliAODHandler") ) { + fAOD = ((AliAODHandler*)handler)->GetAOD(); + } else { + AliFatal("I can't get any AOD Event Handler"); + } + } +} + +//____________________________________________________________________ +void AliAnalysisTaskUE::CreateOutputObjects() +{ +// Create the output container +// + if (fDebug > 1) AliInfo("CreateOutPutData()"); + // + // Histograms + CreateHistos(); + fListOfHistos->SetOwner(kTRUE); +// OpenFile(0); +} + + +//____________________________________________________________________ +void AliAnalysisTaskUE::Exec(Option_t */*option*/) +{ +// Execute analysis for current event +// + AnalyseUE(); + // Post the data + PostData(0, fListOfHistos); +} + + +//____________________________________________________________________ +void AliAnalysisTaskUE::Terminate(Option_t */*option*/) +{ +// Terminate analysis +// + + // Normalize histos to region area TODO: + Double_t areafactor = 1.; + if( fRegionType == 1 ) { + areafactor = 1./ (fTrackEtaCut *2. * TMath::Pi()*2.); + } else { + areafactor = 1./ ( fConeRadius * fConeRadius * TMath::Pi() ); + } + + // Draw some Test plot to the screen + if (!gROOT->IsBatch()) { + + // Update pointers reading them from the output slot + fListOfHistos = dynamic_cast (GetOutputData(0)); + if( !fListOfHistos ) { + AliError("Histogram List is not available"); + return; + } + fhEleadingPt = (TH1F*)fListOfHistos->At(1); + fhdNdEta_PhiDist = (TH1F*)fListOfHistos->At(8); + fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11); + fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12); + fhRegionMultMaxVsEt = (TH1F*)fListOfHistos->At(14); + fhRegionMultMinVsEt = (TH1F*)fListOfHistos->At(15); + fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16); + + fhNJets = (TH1F*)fListOfHistos->At(0); + + // Canvas + TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1200,800); + c1->Divide(2,2); + c1->cd(1); + TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1); + h1r->Scale( areafactor ); + h1r->SetMarkerStyle(20); + h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); + h1r->SetYTitle("P_{T}^{90, max}"); + h1r->DrawCopy("p"); + + c1->cd(2); + + TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1); + h2r->Scale( areafactor ); + h2r->SetMarkerStyle(20); + h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); + h2r->SetYTitle("P_{T}^{90, min}"); + h2r->DrawCopy("p"); + + c1->cd(3); + TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1); + h4r->Scale(2.); // make average + h4r->Scale( areafactor ); + h4r->SetYTitle("#DeltaP_{T}^{90}"); + h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); + h4r->SetMarkerStyle(20); + h4r->DrawCopy("p"); + + c1->cd(4); + TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + TH1F *h6r = new TH1F("hRegionMultMinVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + + h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1); + h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1); + h5r->Scale( areafactor ); + h6r->Scale( areafactor ); + h5r->SetYTitle("N_{Tracks}^{90}"); + h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); + h5r->SetMarkerStyle(20); + h5r->DrawCopy("p"); + h6r->SetMarkerStyle(21); + h6r->SetMarkerColor(2); + h6r->SetYTitle("N_{Tracks}^{90}"); + h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); + h6r->DrawCopy("p same"); + c1->Update(); + + TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800); + c2->Divide(2,2); + c2->cd(1); + fhEleadingPt->SetMarkerStyle(20); + fhEleadingPt->SetMarkerColor(2); + fhEleadingPt->DrawCopy("p"); + gPad->SetLogy(); + + c2->cd(2); + fhdNdEta_PhiDist->SetMarkerStyle(20); + fhdNdEta_PhiDist->SetMarkerColor(2); + fhdNdEta_PhiDist->DrawCopy("p"); + gPad->SetLogy(); + + c2->cd(3); + fhNJets->DrawCopy(); + + //fhTransRegPartPtDist = (TH1F*)fListOfHistos->At(2); + fhRegionMultMin = (TH1F*)fListOfHistos->At(3); + fhMinRegAvePt = (TH1F*)fListOfHistos->At(4); + fhMinRegSumPt = (TH1F*)fListOfHistos->At(5); + //fhMinRegMaxPtPart = (TH1F*)fListOfHistos->At(6); + fhMinRegSumPtvsMult = (TH1F*)fListOfHistos->At(7); + + // Canvas + TCanvas* c3 = new TCanvas("c3"," p_{T} dist",160,160,1200,800); + c3->Divide(2,2); + c3->cd(1); + /*fhTransRegPartPtDist->SetMarkerStyle(20); + fhTransRegPartPtDist->SetMarkerColor(2); + fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries()); + fhTransRegPartPtDist->DrawCopy("p"); + gPad->SetLogy(); + */ + c3->cd(2); + fhMinRegSumPt->SetMarkerStyle(20); + fhMinRegSumPt->SetMarkerColor(2); + fhMinRegSumPt->Scale(areafactor); + fhMinRegSumPt->DrawCopy("p"); + gPad->SetLogy(); + + c3->cd(3); + fhMinRegAvePt->SetMarkerStyle(20); + fhMinRegAvePt->SetMarkerColor(2); + fhMinRegAvePt->Scale(areafactor); + fhMinRegAvePt->DrawCopy("p"); + gPad->SetLogy(); + + c3->cd(4); + + TH1F *h7r = new TH1F("hRegionMultMinVsMult", "", 21, -0.5, 20.5); + + h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1); + + h7r->SetMarkerStyle(20); + h7r->SetMarkerColor(2); + h7r->DrawCopy("p"); + + c3->Update(); + + +/* c2->cd(3); + fhValidRegion->SetMarkerStyle(20); + fhValidRegion->SetMarkerColor(2); + fhValidRegion->DrawCopy("p"); +*/ + c2->Update(); + } else { + AliInfo(" Batch mode, not histograms will shown..."); + } + + if( fDebug > 1 ) + AliInfo("End analysis"); + +} + + +//____________________________________________________________________ +void AliAnalysisTaskUE::AnalyseUE() +{ + // + // Look for correlations on the tranverse regions to + // the leading charged jet + // + // ------------------------------------------------ + // Find Leading Jets 1,2,3 + // (could be skipped if Jets are sort by Pt...) + Double_t maxPtJet1 = 0.; + Int_t index1 = -1; + Double_t maxPtJet2 = 0.; // jet 2 need for back to back inclusive + Int_t index2 = -1; + Double_t maxPtJet3 = 0.; // jet 3 need for back to back exclusive + Int_t index3 = -1; + Int_t nJets = fAOD->GetNJets(); + for( Int_t i=0; iGetJet(i); + if( jet->Pt() > maxPtJet1 ) { + maxPtJet3 = maxPtJet2; index3 = index2; + maxPtJet2 = maxPtJet1; index2 = index1; + maxPtJet1 = jet->Pt(); index1 = i; + } else if( jet->Pt() > maxPtJet2 ) { + maxPtJet3 = maxPtJet2; index3 = index2; + maxPtJet2 = jet->Pt(); index2 = i; + } else if( jet->Pt() > maxPtJet3 ) { + maxPtJet3 = jet->Pt(); index3 = i; + } + } + + if( index1 < 0 ) { + AliInfo("\nSkipping Event, not jet found..."); + return; + } else + AliInfo(Form("\nPt Leading Jet = %6.1f eta=%5.1f ", maxPtJet1, fAOD->GetJet(index1)->Eta() )); + + fhNJets->Fill(nJets); + + TVector3 jetVect[2]; + + // ---------------------------------------------- + // Cut events by jets topology + // fAnaType: + // 1 = inclusive, + // - Jet1 |eta| < fJet1EtaCut + // 2 = back to back inclusive + // - fulfill case 1 + // - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut + // - Jet2.Pt/Jet1Pt > fJet2RatioPtCut + // 3 = back to back exclusive + // - fulfill case 2 + // - Jet3.Pt < fJet3PtCut + + Bool_t isInclusive = kFALSE; + + if( TMath::Abs(fAOD->GetJet(index1)->Eta()) > fJet1EtaCut) { + if( fDebug > 1 ) AliInfo("Skipping Event...Jet1 |eta| > fJet1EtaCut"); + return; + } + isInclusive = kTRUE; + jetVect[0].SetXYZ(fAOD->GetJet(index1)->Px(), + fAOD->GetJet(index1)->Py(), + fAOD->GetJet(index1)->Pz()); + // back to back inclusive + Bool_t isB2Binclusive = kFALSE; + if( fAnaType > 1 && index2 > 0 && isInclusive) { + jetVect[1].SetXYZ(fAOD->GetJet(index2)->Px(), + fAOD->GetJet(index2)->Py(), + fAOD->GetJet(index2)->Pz()); + if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut || + maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) { + if( fDebug > 1 ) AliInfo("Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut"); + return; + } + isB2Binclusive = kTRUE; + } + if (isInclusive && !isB2Binclusive && fAnaType>1) return; + // back to back exclusive + Bool_t isB2Bexclusive = kFALSE; + if( fAnaType > 2 && index3 > 0 && isB2Binclusive) { + if( fAOD->GetJet(index3)->Pt() > fJet3PtCut ) { + if( fDebug > 1 ) AliInfo("Skipping Event... Jet3.Pt > fJet3PtCut "); + return; + } + isB2Bexclusive = kTRUE; + } + if (isB2Binclusive && !isB2Bexclusive && fAnaType>2) return; + + AliInfo(Form("njets = %d",nJets)); + fhEleadingPt->Fill( maxPtJet1 ); + + // ---------------------------------------------- + // Find max and min regions + Double_t sumPtRegionPosit = 0.; + Double_t sumPtRegionNegat = 0.; + Double_t maxPartPtRegion = 0.; + Int_t nTrackRegionPosit = 0; + Int_t nTrackRegionNegat = 0; + static const Double_t k270rad = 270.*TMath::Pi()/180.; + + Int_t nTracks = fAOD->GetNTracks(); + for (Int_t ipart=0; ipartGetTrack( ipart ); + + if ( !part->Charge() ) continue; + if ( part->Pt() < fTrackPtCut ) continue; + + TVector3 partVect(part->Px(), part->Py(), part->Pz()); + + Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad; + if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi(); + fhdNdEta_PhiDist->Fill( deltaPhi ); + fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); + + Int_t region = IsTrackInsideRegion( jetVect, &partVect ); + + if (region > 0) { + if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt(); + sumPtRegionPosit += part->Pt(); + nTrackRegionPosit++; + fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); + } + if (region < 0) { + if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt(); + sumPtRegionNegat += part->Pt(); + nTrackRegionNegat++; + fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); + } + } + + //How quantities will be sorted before Fill Min and Max Histogram + // 1=Plots will be CDF-like + // 2=Plots will be Marchesini-like + if (fOrdering==1){ + if (sumPtRegionPosit > sumPtRegionNegat) { + FillSumPtRegion( maxPtJet1, sumPtRegionPosit, sumPtRegionNegat ); + } else { + FillSumPtRegion( maxPtJet1, sumPtRegionNegat, sumPtRegionPosit ); + } + if (nTrackRegionPosit > nTrackRegionNegat ) { + FillMultRegion( maxPtJet1, nTrackRegionPosit, nTrackRegionNegat, sumPtRegionNegat ); + } else { + FillMultRegion( maxPtJet1, nTrackRegionNegat, nTrackRegionPosit, sumPtRegionPosit ); + } + } else if (fOrdering==2) { + if (sumPtRegionPosit > sumPtRegionNegat) { + FillSumPtRegion( maxPtJet1, sumPtRegionPosit, sumPtRegionNegat ); + FillMultRegion( maxPtJet1, nTrackRegionPosit, nTrackRegionNegat, sumPtRegionNegat ); + } else { + FillSumPtRegion( maxPtJet1, sumPtRegionNegat, sumPtRegionPosit ); + FillMultRegion( maxPtJet1, nTrackRegionNegat, nTrackRegionPosit, sumPtRegionPosit ); + } + } + + Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.; + Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.; + if (avePosRegion > aveNegRegion) { + FillAvePartPtRegion( maxPtJet1, avePosRegion, aveNegRegion ); + } else { + FillAvePartPtRegion( maxPtJet1, aveNegRegion, avePosRegion ); + } + + fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion ); + + // Compute pedestal like magnitude + fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/2.0); + fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/2.0); + +} + + +//____________________________________________________________________ +void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ) +{ + // Fill sumPt of control regions + + // Max cone + fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax ); + // Min cone + fhRegionSumPtMinVsEt->Fill( leadingE, ptMin ); + // MAke distributions for UE comparison with MB data + fhMinRegSumPt->Fill(ptMin); +} + +//____________________________________________________________________ +void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ) +{ + // Fill average particle Pt of control regions + + // Max cone + fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax ); + // Min cone + fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin ); + // MAke distributions for UE comparison with MB data + fhMinRegAvePt->Fill(ptMin); +} + +//____________________________________________________________________ +void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin ) +{ + // Fill Nch multiplicity of control regions + + // Max cone + fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax ); + fhRegionMultMax->Fill( nTrackPtmax ); + // Min cone + fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin ); + fhRegionMultMin->Fill( nTrackPtmin ); + // MAke distributions for UE comparison with MB data + fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin); + +} + +//____________________________________________________________________ +Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect) +{ + // return de region in delta phi + // -1 negative delta phi + // 1 positive delta phi + // 0 outside region + static const Double_t k60rad = 60.*TMath::Pi()/180.; + static const Double_t k120rad = 120.*TMath::Pi()/180.; + + Int_t region = 0; + if( fRegionType == 1 ) { + if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0; + // transverse regions + if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; + if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1; + + } else if( fRegionType == 2 ) { + // Cone regions + Double_t deltaR = 0.; + + TVector3 positVect,negatVect; + positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2()); + negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2()); + + if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) { + region = 1; + deltaR = positVect.DrEtaPhi(*partVect); + } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) { + region = -1; + deltaR = negatVect.DrEtaPhi(*partVect); + } + + if (deltaR > fConeRadius) region = 0; + + } else + AliError("Unknow region type"); + + // For debug (to be removed) + //if( region != 0 ) fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) ); + + return region; +} + +//____________________________________________________________________ +void AliAnalysisTaskUE::CreateHistos() +{ + fListOfHistos = new TList(); + + + fhNJets = new TH1F("fhNJets", "n Jet", 10, 0, 10); + fhNJets->SetXTitle("# of jets"); + fhNJets->Sumw2(); + fhNJets->Sumw2(); + fListOfHistos->Add( fhNJets ); // At(0) + + fhEleadingPt = new TH1F("hEleadingPt", "Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhEleadingPt->SetXTitle("P_{T} (GeV/c)"); + fhEleadingPt->SetYTitle("dN/dP_{T}"); + fhEleadingPt->Sumw2(); + fListOfHistos->Add( fhEleadingPt ); // At(1) + + fhMinRegPtDist = new TH1F("hMinRegPtDist", "P_{T} distribution in Min zone", 50,0.,20.); + fhMinRegPtDist->SetXTitle("P_{T} (GeV/c)"); + fhMinRegPtDist->SetYTitle("dN/dP_{T}"); + fhMinRegPtDist->Sumw2(); + fListOfHistos->Add( fhMinRegPtDist ); // At(2) + + fhRegionMultMin = new TH1F("hRegionMultMin", "N_{ch}^{90, min}", 21, -0.5, 20.5); + fhRegionMultMin->SetXTitle("N_{ch tracks}"); + fhRegionMultMin->Sumw2(); + fListOfHistos->Add( fhRegionMultMin ); // At(3) + + fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT", 50, 0., 20.); + fhMinRegAvePt->SetXTitle("P_{T} (GeV/c)"); + fhMinRegAvePt->Sumw2(); + fListOfHistos->Add( fhMinRegAvePt ); // At(4) + + fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ", 50, 0., 20.); + fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})"); + fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)"); + fhMinRegSumPt->Sumw2(); + fListOfHistos->Add( fhMinRegSumPt ); // At(5) + + fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ", 50, 0., 20.); + fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})"); + fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)"); + fhMinRegMaxPtPart->Sumw2(); + fListOfHistos->Add( fhMinRegMaxPtPart ); // At(6) + + fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ", 21, -0.5, 20.5); + fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)"); + fhMinRegSumPtvsMult->SetXTitle("N_{charge}"); + fhMinRegSumPtvsMult->Sumw2(); + fListOfHistos->Add( fhMinRegSumPtvsMult ); // At(7); + + fhdNdEta_PhiDist = new TH1F("hdNdEta_PhiDist", "Charge particle density |#eta|< 1 vs #Delta#phi", 120, 0., 2.*TMath::Pi()); + fhdNdEta_PhiDist->SetXTitle("#Delta#phi"); + fhdNdEta_PhiDist->SetYTitle("dN_{ch}/d#etad#phi"); + fhdNdEta_PhiDist->Sumw2(); + fListOfHistos->Add( fhdNdEta_PhiDist ); // At(8) + + // Can be use to get part pt distribution for differente Jet Pt bins + fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", "dN/dP_{T} |#eta|<1 vs Leading Jet P_{T}", + 50,0.,50., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}"); + fhFullRegPartPtDistVsEt->SetXTitle("p_{T}"); + fhFullRegPartPtDistVsEt->Sumw2(); + fListOfHistos->Add( fhFullRegPartPtDistVsEt ); // At(9) + + // Can be use to get part pt distribution for differente Jet Pt bins + fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", "dN/dP_{T} in tranvese regions |#eta|<1 vs Leading Jet P_{T}", + 50,0.,50., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}"); + fhTransRegPartPtDistVsEt->SetXTitle("p_{T}"); + fhTransRegPartPtDistVsEt->Sumw2(); + fListOfHistos->Add( fhTransRegPartPtDistVsEt ); // At(10) + + + fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt", "P_{T}^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)"); + fhRegionSumPtMaxVsEt->Sumw2(); + fListOfHistos->Add( fhRegionSumPtMaxVsEt ); // At(11) + + fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt", "P_{T}^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)"); + fhRegionSumPtMinVsEt->Sumw2(); + fListOfHistos->Add( fhRegionSumPtMinVsEt ); // At(12) + + fhRegionMultMax = new TH1I("hRegionMultMax", "N_{ch}^{90, max}", 21, -0.5, 20.5); + fhRegionMultMax->SetXTitle("N_{ch tracks}"); + fhRegionMultMax->Sumw2(); + fListOfHistos->Add( fhRegionMultMax ); // At(13) + + fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt", "N_{ch}^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)"); + fhRegionMultMaxVsEt->Sumw2(); + fListOfHistos->Add( fhRegionMultMaxVsEt ); // At(14) + + fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt", "N_{ch}^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhRegionMultMinVsEt->SetXTitle("E (GeV/c)"); + fhRegionMultMinVsEt->Sumw2(); + fListOfHistos->Add( fhRegionMultMinVsEt ); // At(15) + + fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)"); + fhRegionAveSumPtVsEt->Sumw2(); + fListOfHistos->Add( fhRegionAveSumPtVsEt ); // At(16) + + fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)"); + fhRegionDiffSumPtVsEt->Sumw2(); + fListOfHistos->Add( fhRegionDiffSumPtVsEt ); // At(17) + + fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)"); + fhRegionAvePartPtMaxVsEt->Sumw2(); + fListOfHistos->Add( fhRegionAvePartPtMaxVsEt ); // At(18) + + fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)"); + fhRegionAvePartPtMinVsEt->Sumw2(); + fListOfHistos->Add( fhRegionAvePartPtMinVsEt ); // At(19) + + fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)"); + fhRegionMaxPartPtMaxVsEt->Sumw2(); + fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt ); // At(20) + +/* + // For debug region selection + fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi", + 20, -2.,2., 62, -TMath::Pi(), TMath::Pi()); + fhValidRegion->SetYTitle("#Delta#phi"); + fhValidRegion->Sumw2(); + fListOfHistos->Add( fhValidRegion ); // At(15) +*/ +} + + + diff --git a/PWG4/AliAnalysisTaskUE.h b/PWG4/AliAnalysisTaskUE.h new file mode 100644 index 00000000000..974fab08ace --- /dev/null +++ b/PWG4/AliAnalysisTaskUE.h @@ -0,0 +1,138 @@ +#ifndef ALIANALYSISTASKUE_H +#define ALIANALYSISTASKUE_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +#include "AliAnalysisTaskSE.h" + +class AliESDEvent; +class AliAODEvent; +class TH1F; +class TH2F; +class TH1I; +class TVector3; + +class AliAnalysisTaskUE : public AliAnalysisTask +{ + public: + AliAnalysisTaskUE(const char* name="AliAnalysisTaskUE"); + virtual ~AliAnalysisTaskUE() {;} + + // Implementation of interface methods + virtual void ConnectInputData(Option_t *); + virtual void CreateOutputObjects(); + virtual void Exec(Option_t *option); + virtual void Terminate(Option_t *); + + // Setters + virtual void SetDebugLevel(Int_t level) { fDebug = level; } + void SetPtRangeInHist(Int_t bin, Double_t min, Double_t max) { + fBinsPtInHist = bin; + fMinJetPtInHist = min; + fMaxJetPtInHist = max; + } + + void SetAnaTopology(Int_t val) { fAnaType = val; } + void SetRegionType(Int_t val) { fRegionType = val; } + void SetConeRadius(Double_t val) { fConeRadius = val; } + void SetPtSumOrdering(Int_t val) { fOrdering = val; } + // Jet cuts + void SetJet1EtaCut(Double_t val) { fJet1EtaCut = val; } + void SetJet2DeltaPhiCut(Double_t val) { fJet2DeltaPhiCut = val; } + void SetJet2RatioPtCut(Double_t val) { fJet2RatioPtCut = val; } + void SetJet3PtCut(Double_t val) { fJet3PtCut = val; } + // track cuts + void SetTrackPtCut(Double_t val) { fTrackPtCut = val; } + void SetTrackEtaCut(Double_t val) { fTrackEtaCut = val; } + +private: + AliAnalysisTaskUE(const AliAnalysisTaskUE &det); +AliAnalysisTaskUE& operator=(const AliAnalysisTaskUE &det); + + void AnalyseUE(); + Int_t IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect); + void CreateHistos(); + void FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ); + void FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ); + void FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin ); + + + Int_t fDebug; // Debug flag + AliAODEvent* fAOD; //! AOD Event + TList* fListOfHistos; // Output list of histograms + + // Config + Int_t fBinsPtInHist; // # bins for Pt histos range + Double_t fMinJetPtInHist; // min Jet Pt value for histo range + Double_t fMaxJetPtInHist; // max Jet Pt value for histo range + + // Cuts + Int_t fAnaType; // Analysis type on jet topology: + // 1=inclusive (default) + // 2=back to back inclusive + // 3=back to back exclusive + // 4=gama jet (back to back) ??? + // Minimum bias + // 31 = Semi jet (charged leading particle jets) + // 32 = Random jetcone ? + // 33 = Swiss chees ? + + Int_t fRegionType; // 1 = transverse regions (default) + // 2 = cone regions + + Double_t fConeRadius; // if selected Cone-like region type set Radius (=0.7 default) + + Int_t fOrdering; // Pt and multiplicity summation ordering: + // 1=CDF-like -independent sorting according quantity to be scored: Double sorting- (default) + // if Pt summation will be scored take Pt minimum between both zones and + // fill Pt Max. and Min. histog. accordingly + // if Multiplicity summation will be scored take Mult. minimum between both zones and + // fill Mult Max and Min histog. accordingly + // 2=Marchesini-like (Only Pt sorting: Single sorting) + // sort only according Pt summation scored, find minimum between both zones and + // fill Pt and Multiplicity Max and Min summation histog. following only this criterium + // 3=User Selection sorting (NOTE: USER must implement it within cxx) + + // Jet cuts + Double_t fJet1EtaCut; // |jet1 eta| < fJet1EtaCut (fAnaType = 1,2,3) + Double_t fJet2DeltaPhiCut; // |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut (fAnaType = 2,3) + Double_t fJet2RatioPtCut; // Jet2.Pt/Jet1Pt > fJet2RatioPtCut (fAnaType = 2,3) + Double_t fJet3PtCut; // Jet3.Pt < fJet3PtCut (fAnaType = 3) + // track cuts + Double_t fTrackPtCut; // Pt cut of tracks in the regions + Double_t fTrackEtaCut; // Eta cut on tracks in the regions (fRegionType=1) + + // Histograms ( are owner by fListOfHistos TList ) + TH1F* fhNJets; //! + TH1F* fhEleadingPt; //! + + TH1F* fhMinRegPtDist; + TH1F* fhRegionMultMin; + TH1F* fhMinRegAvePt; + TH1F* fhMinRegSumPt; + TH1F* fhMinRegMaxPtPart; + TH1F* fhMinRegSumPtvsMult; + + TH1F* fhdNdEta_PhiDist; //! + TH2F* fhFullRegPartPtDistVsEt; //! + TH2F* fhTransRegPartPtDistVsEt; //! + + TH1F* fhRegionSumPtMaxVsEt; //! + TH1I* fhRegionMultMax; //! + TH1F* fhRegionMultMaxVsEt; //! + TH1F* fhRegionSumPtMinVsEt; //! + TH1F* fhRegionMultMinVsEt; //! + TH1F* fhRegionAveSumPtVsEt; //! + TH1F* fhRegionDiffSumPtVsEt; //! + + TH1F* fhRegionAvePartPtMaxVsEt; //! + TH1F* fhRegionAvePartPtMinVsEt; //! + TH1F* fhRegionMaxPartPtMaxVsEt; //! + + // TH2F* fhValidRegion; //! test to be canceled + + ClassDef( AliAnalysisTaskUE, 1); // Analysis task for Underlying Event analysis +}; + +#endif diff --git a/PWG4/Makefile b/PWG4/Makefile index 2ab69a08578..ccd2932d0d4 100644 --- a/PWG4/Makefile +++ b/PWG4/Makefile @@ -1,38 +1,58 @@ +PACKAGE = invalid-only-for-proof -include Makefile.arch +include $(ROOTSYS)/test/Makefile.arch +include lib$(PACKAGE).pkg -default-target: libPWG4PartCorr.so +ifndef PACKCXXFLAGS + PACKCXXFLAGS = $(CXXFLAGS) +endif -ALICEINC = -I. +ALICEINC = -I. ### define include dir for local case and par case +ifneq ($(STEERBase_INCLUDE),) + ALICEINC += -I../$(STEERBase_INCLUDE) +endif + +ifneq ($(ESD_INCLUDE),) + ALICEINC += -I../$(ESD_INCLUDE) +endif + +ifneq ($(AOD_INCLUDE),) + ALICEINC += -I../$(AOD_INCLUDE) +endif + ifneq ($(ANALYSIS_INCLUDE),) - ALICEINC += -I../$(ESD_INCLUDE) -I../$(AOD_INCLUDE) -I../$(ANALYSIS_INCLUDE) -I../$(STEERBase_INCLUDE) -I../$(ANALYSISalice_INCLUDE) -else - ifneq ($(ALICE_ROOT),) - ALICEINC += -I$(ALICE_ROOT)/include - endif + ALICEINC += -I../$(ANALYSIS_INCLUDE) endif -# for building of PWG4PartCorr.par -ifneq ($(GAMMA_INCLUDE),) - ALICEINC += -I../$(GAMMA_INCLUDE) +ifneq ($(ANALYSISalice_INCLUDE),) + ALICEINC += -I../$(ANALYSISalice_INCLUDE) endif -CXXFLAGS += $(ALICEINC) -g +ifneq ($(PWG4PartCorr_INCLUDE),) + ALICEINC += -I../$(PWG4PartCorr_INCLUDE) +endif + + +ifneq ($(PWG4JetTasks_INCLUDE),) + ALICEINC += -I../$(PWG4JetTasks_INCLUDE) +endif -PACKAGE = PWG4PartCorr -include lib$(PACKAGE).pkg -DHDR_PWG4PartCorr := $(DHDR) -HDRS_PWG4PartCorr := $(HDRS) -SRCS_PWG4PartCorr := $(SRCS) G__$(PACKAGE).cxx -OBJS_PWG4PartCorr := $(SRCS_PWG4PartCorr:.cxx=.o) +# only if no par file was loaded before +ifeq ($(ALICEINC),-I.) + ifneq ($(ALICE_ROOT),) + ALICEINC += -I$(ALICE_ROOT)/include + endif +endif -PARFILE = $(PACKAGE).par +CXXFLAGS += $(ALICEINC) -g +SRCS += G__$(PACKAGE).cxx +OBJS = $(SRCS:.cxx=.o) -lib$(PACKAGE).so: $(OBJS_PWG4PartCorr) +lib$(PACKAGE).so: $(OBJS) @echo "Linking" $@ ... @/bin/rm -f $@ ifeq ($(ARCH),macosx) @@ -44,46 +64,12 @@ endif @echo "done" %.o: %.cxx %.h - $(CXX) $(CXXFLAGS) -c $< -o $@ + $(CXX) $(PACKCXXFLAGS) -c $< -o $@ clean: - @rm -f $(OBJS_PWG4PartCorr) *.so G__$(PACKAGE).* $(PARFILE) + @rm -f $(OBJS) *.so G__$(PACKAGE).* G__$(PACKAGE).cxx G__$(PACKAGE).h: $(HDRS) $(DHDR) - @echo "Generating dictionary ..." - rootcint -f $@ -c $(ALICEINC) $^ - -### CREATE PAR FILE - -$(PARFILE): $(patsubst %,$(PACKAGE)/%,$(filter-out G__%, $(HDRS_PWG4PartCorr) $(SRCS_PWG4PartCorr) $(DHDR_PWG4PartCorr) Makefile Makefile.arch lib$(PACKAGE).pkg PROOF-INF)) - @echo "Creating archive" $@ ... - @tar cfzh $@ $(PACKAGE) - @rm -rf $(PACKAGE) - @echo "done" + @echo "Generating dictionaries ..." $(ALICEINC) + rootcint -f $@ -c $(CINTFLAGS) $(ALICEINC) $^ -$(PACKAGE)/Makefile: Makefile #.$(PACKAGE) - @echo Copying $< to $@ with transformations - @[ -d $(dir $@) ] || mkdir -p $(dir $@) - @sed 's/include \$$(ROOTSYS)\/test\/Makefile.arch/include Makefile.arch/' < $^ > $@ - -$(PACKAGE)/Makefile.arch: $(ROOTSYS)/test/Makefile.arch - @echo Copying $< to $@ - @[ -d $(dir $@) ] || mkdir -p $(dir $@) - @cp -a $^ $@ - -$(PACKAGE)/PROOF-INF: PROOF-INF.$(PACKAGE) - @echo Copying $< to $@ - @[ -d $(dir $@) ] || mkdir -p $(dir $@) - @cp -a -r $^ $@ - -$(PACKAGE)/%: % - @echo Copying $< to $@ - @[ -d $(dir $@) ] || mkdir -p $(dir $@) - @cp -a $< $@ - -test-%.par: %.par - @echo "INFO: The file $< is now tested, in case of an error check in par-tmp." - @mkdir -p par-tmp - @cd par-tmp; tar xfz ../$<; cd $(subst .par,,$<); PROOF-INF/BUILD.sh - @rm -rf par-tmp - @echo "INFO: Testing succeeded (already cleaned up)" diff --git a/PWG4/PROOF-INF.PWG4JetTasks/BUILD.sh b/PWG4/PROOF-INF.PWG4JetTasks/BUILD.sh new file mode 100755 index 00000000000..fc9490a6c2d --- /dev/null +++ b/PWG4/PROOF-INF.PWG4JetTasks/BUILD.sh @@ -0,0 +1,3 @@ +#! /bin/sh + +make diff --git a/PWG4/PROOF-INF.PWG4JetTasks/SETUP.C b/PWG4/PROOF-INF.PWG4JetTasks/SETUP.C new file mode 100644 index 00000000000..1d79478be01 --- /dev/null +++ b/PWG4/PROOF-INF.PWG4JetTasks/SETUP.C @@ -0,0 +1,13 @@ +void SETUP() +{ + + // Load the JET-Tasks library + gSystem->Load("libPWG4JetTasks"); + + // Set the Include paths + gSystem->AddIncludePath("-I$ROOTSYS/include -IPWG4"); + gROOT->ProcessLine(".include PWG4"); + + // Set our location, so that other packages can find us + gSystem->Setenv("PWG4JetTasks_INCLUDE", "PWG4"); +} diff --git a/PWG4/PWG4JetTasksLinkDef.h b/PWG4/PWG4JetTasksLinkDef.h new file mode 100644 index 00000000000..80b9442efd3 --- /dev/null +++ b/PWG4/PWG4JetTasksLinkDef.h @@ -0,0 +1,10 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class AliAnalysisTaskUE+; + + +#endif diff --git a/PWG4/libPWG4JetTasks.pkg b/PWG4/libPWG4JetTasks.pkg new file mode 100644 index 00000000000..b41ddad5d07 --- /dev/null +++ b/PWG4/libPWG4JetTasks.pkg @@ -0,0 +1,15 @@ +#-*- Mode: Makefile -*- + +SRCS = AliAnalysisTaskUE.cxx + +HDRS:= $(SRCS:.cxx=.h) + +DHDR:= PWG4JetTasksLinkDef.h + +EXPORT:=$(SRCS:.cxx=.h) + +ifeq (win32gcc,$(ALICE_TARGET)) +PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) -lSTEERBase \ + -lESD -lAOD -lANALYSIS -lANALYSISalice \ + -L$(shell root-config --libdir) -lEG +endif diff --git a/PWG4/macros/AnalysisTrainCAF.C b/PWG4/macros/AnalysisTrainCAF.C new file mode 100644 index 00000000000..27f62bdf7c8 --- /dev/null +++ b/PWG4/macros/AnalysisTrainCAF.C @@ -0,0 +1,256 @@ +//______________________________________________________________________________ +void AnalysisTrainCAF() +{ +// Example of running analysis train in CAF. To run in debug mode: +// - export ROOTSYS=debug on your local client +// - un-comment gProof->ClearPackages() +// - un-comme lines with debugging info + +// WHY AOD is not a exchange container when running from ESD->AOD + + Bool_t debug = kTRUE; + Bool_t useMC = kTRUE; + Bool_t readTR = kFALSE; + + Int_t iAODanalysis = 0; + Int_t iAODhandler = 1; + Int_t iESDfilter = 1; // Only active if iAODanalysis=0 + Int_t iJETAN = 1; + Int_t iJETANMC = 1; + Int_t iPWG4UE = 1; + + if (iAODanalysis) { + useMC = kFALSE; + readTR = kFALSE; + iESDfilter = 0; + iMUONfilter = 0; + } + if (iJETAN) iESDfilter=1; + if (iESDfilter) iAODhandler=1; + + // Dataset from CAF + TString dataset = "/PWG4/arian/jetjet15-50"; + printf("==================================================================\n"); + printf("=========== RUNNING ANALYSIS TRAIN IN CAF MODE =============\n"); + printf("==================================================================\n"); + if (iAODanalysis) printf("= AOD analysis on dataset: %s\n", dataset.Data()); + else printf("= ESD analysis on dataset: %s\n", dataset.Data()); + if (iESDfilter) printf("= ESD filter =\n"); + if (iJETAN) printf("= Jet analysis =\n"); + if (iJETANMC) printf("= Jet analysis from Kinematics =\n"); + if (iPWG4UE) printf("= PWG4 UE =\n"); + printf("==================================================================\n"); + if (useMC) printf(":: use MC TRUE\n"); + else printf(":: use MC FALSE\n"); + if (readTR) printf(":: read TR TRUE\n"); + else printf(":: read TR FALSE\n"); + if (debug) printf(":: debugging TRUE\n"); + else printf(":: debugging FALSE\n"); + + // Load common libraries + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libPhysics.so"); + + + // Reset user processes if CAF if not responding anymore + //TProof::Reset("lxb6046"); + // One may enable a different ROOT version on CAF + + const char* proofNode = "lxb6046"; + // const char* proofNode = "kleinb@localhost"; + + TProof::Mgr(proofNode)->ShowROOTVersions(); + // TProof::Mgr(proofNode)->SetROOTVersion("vHEAD-r24503_dbg"); + + // Connect to proof + TProof::Open(proofNode); // may be username@lxb6046 if user not the same as on local + + // Clear packages if changing ROOT version on CAF or local + // gProof->ClearPackages(); + // Enable proof debugging if needed + // gProof->SetLogLevel(5); + // To debug the train in PROOF mode, type in a root session: + // root[0] TProof::Mgr("lxb6064")->GetSessionLogs()->Display("*",0,10000); + // Common packages + // --- Enable the STEERBase Package + gProof->UploadPackage("pars/STEERBase.par"); + gProof->EnablePackage("STEERBase"); + // --- Enable the ESD Package + gProof->UploadPackage("pars/ESD.par"); + gProof->EnablePackage("ESD"); + // --- Enable the AOD Package + gProof->UploadPackage("pars/AOD.par"); + gProof->EnablePackage("AOD"); + // --- Enable the ANALYSIS Package + gProof->UploadPackage("pars/ANALYSIS.par"); + gProof->EnablePackage("ANALYSIS"); + // --- Enable the ANALYSISalice Package + gProof->UploadPackage("pars/ANALYSISalice.par"); + gProof->EnablePackage("ANALYSISalice"); + + AliPDG::AddParticlesToPdgDataBase(); + + // --- Enable the JETAN Package + if (iJETAN||iJETANMC) { + gProof->UploadPackage("pars/JETAN.par"); + gProof->EnablePackage("JETAN"); + } + // --- Enable particle correlation analysis + if (iPWG4UE) { + gProof->UploadPackage("pars/PWG4JetTasks.par"); + gProof->EnablePackage("PWG4JetTasks"); + } + + // Create the chain + // + + // Make the analysis manager + AliAnalysisManager *mgr = new AliAnalysisManager("Analysis Train", "A test setup for the analysis train"); + // Top container for input + AliAnalysisDataContainer *cinput = mgr->CreateContainer("cInput",TChain::Class(), + AliAnalysisManager::kInputContainer); + if (iAODanalysis) { + // AOD input handler + AliAODInputHandler *aodH = new AliAODInputHandler(); + mgr->SetInputEventHandler(aodH); + } else { + // ESD input handler + AliESDInputHandler *esdHandler = new AliESDInputHandler(); + mgr->SetInputEventHandler(esdHandler); +// esdHandler->SetInactiveBranches("FMD CaloCluster"); + } + // Monte Carlo handler + if (useMC && !iAODanalysis) { + AliMCEventHandler* mcHandler = new AliMCEventHandler(); + mgr->SetMCtruthEventHandler(mcHandler); + mcHandler->SetReadTR(readTR); + } + + // This container is managed by the AOD handler + AliAnalysisDataContainer *cout_aod = 0; + if (iAODhandler) { + // AOD output handler + AliAODHandler* aodHandler = new AliAODHandler(); + mgr->SetOutputEventHandler(aodHandler); + aodHandler->SetOutputFileName("AliAODs.root"); +// aodHandler->SetCreateNonStandardAOD(); + cout_aod = mgr->CreateContainer("cAOD", TTree::Class(), + AliAnalysisManager::kOutputContainer, "default"); + cout_aod->SetSpecialOutput(); + } + + // Debugging if needed + if (debug) mgr->SetDebugLevel(3); +// AliLog::EnableDebug(kTRUE); +// AliLog::SetGlobalLogLevel(2); + + + if (iESDfilter && !iAODanalysis) { + // Set of cuts plugged into the ESD filter + // + // standard + AliESDtrackCuts* esdTrackCutsL = new AliESDtrackCuts("AliESDtrackCuts", "Loose"); + esdTrackCutsL->SetMinNClustersTPC(50); + esdTrackCutsL->SetMaxChi2PerClusterTPC(3.5); + esdTrackCutsL->SetMaxCovDiagonalElements(2,2,0.5,0.5,2); + esdTrackCutsL->SetRequireTPCRefit(kTRUE); + esdTrackCutsL->SetMinNsigmaToVertex(3); + esdTrackCutsL->SetRequireSigmaToVertex(kTRUE); + esdTrackCutsL->SetAcceptKingDaughters(kFALSE); + // + AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter"); + trackFilter->AddCuts(esdTrackCutsL); + // + // ESD filter task putting standard info to output AOD (with cuts) + AliAnalysisTaskESDfilter *esdfilter = new AliAnalysisTaskESDfilter("ESD Filter"); + esdfilter->SetTrackFilter(trackFilter); + esdfilter->SetDebugLevel(10); + mgr->AddTask(esdfilter); + // Connect to data containers + mgr->ConnectInput (esdfilter, 0, cinput ); + mgr->ConnectOutput (esdfilter, 0, cout_aod ); + } + // Jet analysis + AliAnalysisDataContainer *c_aodjet = 0; + if (iJETAN && !iAODanalysis) { + AliAnalysisTaskJets *jetana = new AliAnalysisTaskJets("JetAnalysis"); + jetana->SetDebugLevel(10); + jetana->SetConfigFile("${ALICE_ROOT}/JETAN/ConfigJetAnalysis.C"); + mgr->AddTask(jetana); + // Output histograms list for jet analysis + AliAnalysisDataContainer *cout_jet = mgr->CreateContainer("jethist", TList::Class(), + AliAnalysisManager::kOutputContainer, "jethist.root"); + // Dummy AOD output container for jet analysis (no client yet) + c_aodjet = mgr->CreateContainer("cAODjet", TTree::Class(), + AliAnalysisManager::kExchangeContainer); + // Connect to data containers + mgr->ConnectInput (jetana, 0, cinput ); + mgr->ConnectOutput (jetana, 0, c_aodjet ); + // mgr->ConnectOutput (jetana, 0, cout_aod ); + mgr->ConnectOutput (jetana, 1, cout_jet ); + } + // Jet analysisMC + AliAnalysisDataContainer *c_aodjetMC = 0; + if (iJETANMC && useMC) { + AliAnalysisTaskJets *jetanaMC = new AliAnalysisTaskJets("JetAnalysisMC"); + jetanaMC->SetDebugLevel(10); + jetanaMC->SetConfigFile("${ALICE_ROOT}/JETAN/ConfigJetAnalysisMC.C"); + jetanaMC->SetNonStdBranch("jetsMC"); + mgr->AddTask(jetanaMC); + // Output histograms list for jet analysis + AliAnalysisDataContainer *cout_jetMC = mgr->CreateContainer("jethistMC", TList::Class(), + AliAnalysisManager::kOutputContainer, "jethistMC.root"); + // Dummy AOD output container for jet analysis (no client yet) + c_aodjetMC = mgr->CreateContainer("cAODjetMC", TTree::Class(), + AliAnalysisManager::kExchangeContainer); + // Connect to data containers + mgr->ConnectInput (jetanaMC, 0, cinput ); + mgr->ConnectOutput (jetanaMC, 0, c_aodjetMC ); + // mgr->ConnectOutput (jetanaMC, 0, cout_aod ); + mgr->ConnectOutput (jetanaMC, 1, cout_jetMC ); + } + + // Particle correlation analysis + if (iPWG4UE) { + AliAnalysisTaskUE* ueana = new AliAnalysisTaskUE("Undelying Event"); + + + // default parameters use a switch via iPWGUE + // or a config file + Int_t anaType =1; + Int_t regType =1; + Double_t jetEtaCut=0.2; + Double_t trackPtCut=0.5; + Double_t trackEtaCut= 0.9; + Double_t rad=0.7; + Double_t deltaPhiCut = 2.616; + + ueana->SetDebugLevel(10); + ueana->SetPtRangeInHist(25, 0., 250.); + ueana->SetAnaTopology(anaType); + ueana->SetRegionType(regType); + ueana->SetJet1EtaCut(jetEtaCut); + ueana->SetTrackPtCut(trackPtCut); + ueana->SetPtSumOrdering(2); + ueana->SetConeRadius(rad); + ueana->SetTrackEtaCut(trackEtaCut); + ueana->SetJet2DeltaPhiCut(deltaPhiCut); + mgr->AddTask(ueana); + + + AliAnalysisDataContainer *coutput1_UE = mgr->CreateContainer("histosUE", TList::Class(),AliAnalysisManager::kOutputContainer, "histosUE.root"); + + // mgr->ConnectInput (ueana, 0, cinput); + mgr->ConnectInput (ueana, 0, c_aodjet); + mgr->ConnectOutput (ueana, 0, coutput1_UE ); + } + + // Run the analysis + // + if (mgr->InitAnalysis()) { + mgr->PrintStatus(); + mgr->StartAnalysis("proof",dataset.Data(), 20000); + } +}