////////////////////////////////////////////////////////////////////////
//
// Analysis class to Correct Underlying Event studies
-// (requires input AOD)
//
-// Different analysis are performed according to input.
+// This class needs as input ESDs.\r
+// The output is an analysis-specific container.\r
//
-// Run on MC:
-// - fraction of diffractive events after different cuts
-// - tracking efficiency and contamination
-//
-// Run on DATA:
-// - fraction of events after different cuts
-// - vertex reconstruction efficiency
-//
-// vallero@physi.uni-heidelberg.de
+// The class is used to get the contamination from secondaries\r
+// from tracks DCA distribution \r
+// as function of track pT and pseudo-rapidity.\r
+// It provides additional information for the corrections \r
+// that can not be retrieved by the AliAnalysisTaskLeadingTackUE\r
+// task, which is running on AODs.\r
//
////////////////////////////////////////////////////////////////////////
#include <TROOT.h>
#include <TChain.h>
-#include <TCanvas.h>
-#include <TFile.h>
+//#include <TCanvas.h>\r
+//#include <TFile.h>\r
#include <TList.h>
#include <TMath.h>
-#include <TProfile.h>
+//#include <TProfile.h>\r
#include <TTree.h>
-#include <TVector3.h>
+//#include <TVector3.h>\r
+#include <TH3F.h>\r
-#include "AliAnalyseUE.h"
+#include "AliAnalyseLeadingTrackUE.h"\r
+#include "AliAnalysisFilter.h"\r
#include "AliAnalysisHelperJetTasks.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTaskCorrectionsUE.h"
#include "AliAODHandler.h"
#include "AliESDHandler.h"
#include "AliESDtrack.h"
+#include "AliESDtrackCuts.h"\r
#include "AliESDEvent.h"
#include "AliAODInputHandler.h"
#include "AliESDInputHandler.h"
#include "AliCFManager.h"
#include "AliGenDPMjetEventHeader.h"
#include "AliGenPythiaEventHeader.h"
-#include "AliHistogramsUE.h"
#include "AliInputEventHandler.h"
#include "AliLog.h"
#include "AliMCEventHandler.h"
#include "AliAnalysisHelperJetTasks.h"
+class TCanvas;\r
+class TFile;\r
+class TProfile;\r
+class TVector3; \r
+\r
ClassImp( AliAnalysisTaskCorrectionsUE)
// Define global pointer
//____________________________________________________________________
AliAnalysisTaskCorrectionsUE:: AliAnalysisTaskCorrectionsUE(const char* name):
AliAnalysisTask(name,""),
-fAnaUE(0x0),
-fAOD(0x0),
+fAnalyseUE(0x0),\r
+fDebug(0),\r
+fESDEvent(0x0),\r
fESDHandler(0x0),
-fAODBranch("jets"),
-fCFManager(0x0),
-fDebug(0),
-fHistosUE(0x0),
fInputHandler(0x0),
fListOfHistos(0x0),
+fMcEvent(0x0),\r
fMcHandler(0x0),
-fMcEvent(0x0),
-fBinsPtInHist(30),
-fIsNorm2Area(kTRUE),
-fMaxJetPtInHist(300.),
-fMinJetPtInHist(0.),
-fConstrainDistance(kTRUE),
-fMinDistance(0.2),
-fSimulateChJetPt(kFALSE),
-fUseAliStack(kTRUE),
-fUseMCParticleBranch(kFALSE),
-fnTracksVertex(3), // QA tracks pointing to principal vertex (= 3 default)
-fZVertex(5.),
-fAnaType(1),
-fConePosition(1),
-fConeRadius(0.7),
-fFilterBit(0xFF),
-fJetsOnFly(kFALSE),
-fRegionType(1),
-fUseChargeHadrons(kFALSE),
-fUseChPartJet(kFALSE),
-fUsePositiveCharge(kTRUE),
-fUseSingleCharge(kFALSE),
-fOrdering(1),
-fChJetPtMin(5.0),
-fJet1EtaCut(0.2),
-fJet2DeltaPhiCut(2.616), // 150 degrees
-fJet2RatioPtCut(0.8),
-fJet3PtCut(15.),
+fMode(0),\r
+fOutCFcont(0x0),\r
+fhEntries(0x0),\r
+fhFakes(0x0),\r
+fhPtMCAll(0x0),\r
+fhPtMCPrim(0x0),\r
+fhPtMCSec(0x0),\r
+fhPtMCPrimFake(0x0),\r
+fhPtMCSecFake(0x0),\r
+fnTracksVertex(1), // QA tracks pointing to principal vertex \r
+fZVertex(10.),\r
+fhVertexContributors(0x0),\r
+fhVertexReso(0x0),\r
fTrackEtaCut(0.9),
fTrackPtCut(0.),
-fAvgTrials(1)
+fEsdTrackCuts(0x0),\r
+fEsdTrackCutsSPD(0x0),\r
+fEsdTrackCutsSDD(0x0),\r
+fEsdTrackCutsDCA(0x0)\r
{
// Default constructor
// Define input and output slots here
DefineInput(0, TChain::Class());
// Output slot #0 writes into a TList container
DefineOutput(0, TList::Class());
- DefineOutput(1, AliCFContainer::Class());
-
-}
-//____________________________________________________________________
-AliAnalysisTaskCorrectionsUE:: AliAnalysisTaskCorrectionsUE(const AliAnalysisTaskCorrectionsUE & original):
-AliAnalysisTask(),
-fAnaUE(original.fAnaUE),
-fAOD(original.fAOD),
-fESDHandler(original.fESDHandler),
-fAODBranch(original.fAODBranch),
-fCFManager(original.fCFManager),
-fDebug(original.fDebug),
-fHistosUE(original.fHistosUE),
-fInputHandler(original.fInputHandler),
-fListOfHistos(original.fListOfHistos),
-fMcHandler(original.fMcHandler),
-fMcEvent(original.fMcEvent),
-fBinsPtInHist(original.fBinsPtInHist),
-fIsNorm2Area(original.fIsNorm2Area),
-fMaxJetPtInHist(original.fMaxJetPtInHist),
-fMinJetPtInHist(original.fMinJetPtInHist),
-fConstrainDistance(original.fConstrainDistance),
-fMinDistance(original.fMinDistance),
-fSimulateChJetPt(original.fSimulateChJetPt),
-fUseAliStack(original.fUseAliStack),
-fUseMCParticleBranch(original.fUseMCParticleBranch),
-fnTracksVertex(original.fnTracksVertex), // QA tracks pointing to principal vertex (= 3 default)
-fZVertex(original.fZVertex),
-fAnaType(original.fAnaType),
-fConePosition(original.fConePosition),
-fConeRadius(original.fConeRadius),
-fFilterBit(original.fFilterBit),
-fJetsOnFly(original.fJetsOnFly),
-fRegionType(original.fRegionType),
-fUseChargeHadrons(original.fUseChargeHadrons),
-fUseChPartJet(original.fUseChPartJet),
-fUsePositiveCharge(original.fUsePositiveCharge),
-fUseSingleCharge(original.fUseSingleCharge),
-fOrdering(original.fOrdering),
-fChJetPtMin(original.fChJetPtMin),
-fJet1EtaCut(original.fJet1EtaCut),
-fJet2DeltaPhiCut(original.fJet2DeltaPhiCut), // 150 degrees
-fJet2RatioPtCut(original.fJet2RatioPtCut),
-fJet3PtCut(original.fJet3PtCut),
-fTrackEtaCut(original.fTrackEtaCut),
-fTrackPtCut(original.fTrackPtCut),
-fAvgTrials(original.fAvgTrials)
-{
- // Copy constructor
}
-
//______________________________________________________________
AliAnalysisTaskCorrectionsUE & AliAnalysisTaskCorrectionsUE::operator = (const AliAnalysisTaskCorrectionsUE & /*source*/)
{
return *this;
}
+\r
+/************** INTERFACE METHODS *****************************/\r
+\r
//______________________________________________________________
Bool_t AliAnalysisTaskCorrectionsUE::Notify()
{
- //
- // Implemented Notify() to read the cross sections
- // and number of trials from pyxsec.root
- // Copy from AliAnalysisTaskJFSystematics
- fAvgTrials = 1;
- TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
- Float_t xsection = 0;
- Float_t trials = 1;
- if(tree){
- TFile *curfile = tree->GetCurrentFile();
- if (!curfile) {
- Error("Notify","No current file");
- return kFALSE;
- }
-
- AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
- fHistosUE->GetXsec()->Fill("<#sigma>",xsection);
-
- // construct average trials
- Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
- if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
- }
return kTRUE;
void AliAnalysisTaskCorrectionsUE::ConnectInputData(Option_t* /*option*/)
{
// Connect the input data
-
- // We need AODs with tracks and jets.
- // Since AODs can either be connected to the InputEventHandler
- // or to the OutputEventHandler ( the AOD is created by a previous task in the train )
- // we need to check where it is and get the pointer to AODEvent in the right way
-
- // Delta AODs are also accepted
-
if (fDebug > 1) AliInfo("ConnectInputData() ");
//Get the input handler
- TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
-
- if( handler && handler->InheritsFrom("AliAODInputHandler") ) { // input AOD
- fAOD = ((AliAODInputHandler*)handler)->GetEvent();
- if(!fJetsOnFly){
- if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
- }else{
- if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler / Jets on-the-fly");
- }
- } else { //output AOD
- handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
- if( handler && handler->InheritsFrom("AliAODHandler") ) {
- fAOD = ((AliAODHandler*)handler)->GetAOD();
- if (!fJetsOnFly){
- if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
- } else {
- if (fDebug > 1) AliInfo(" ==== Tracks from AliAODHandler / Jets on-the-fly");
- }
- }else {
- AliFatal("I can't get any AOD Event Handler");
- return;
- }
- }
-
- //Get ESD input handler
- AliInputEventHandler *inputHandler = dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
- AliESDEvent *esdEvent = (AliESDEvent*)inputHandler->GetEvent();
- if (!esdEvent && fDebug > 1) {
- AliInfo("********************** No ESD event: cannot retrive DCA values from AOD !!! ");
- }else fAnaUE->SetESDEvent(esdEvent);
+ fESDHandler = (AliESDInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
+ if ( !fESDHandler && fDebug > 0 ) {\r
+ AliFatal(" No ESD event handler connected !!! ");\r
+ return;\r
+ }\r
+
+ //Get ESD event\r
+ fESDEvent = (AliESDEvent*)fESDHandler->GetEvent();\r
+ if (!fESDEvent && fDebug > 1) {\r
+ AliFatal(" No ESD event retrieved !!! ");\r
+ return;\r
+ }\r
//Get MC handler
fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-
- //Initialize AliAnalysisUE class
- fAnaUE->Initialize(fAnaType, fAOD, fConeRadius, fDebug, fFilterBit, fJet1EtaCut, fJet2DeltaPhiCut, fJet2RatioPtCut, fJet3PtCut, fOrdering, fRegionType, fSimulateChJetPt, fTrackEtaCut, fTrackPtCut, fUseChargeHadrons, fUseChPartJet, fUsePositiveCharge, fUseSingleCharge, fHistosUE);
-
+ \r
+ // Define track cuts\r
+ fEsdTrackCuts = new AliESDtrackCuts("ITSTPC", "ITS+TPC standard 2009 cuts w.o. SPD cluster requirement nor DCA cut");\r
+ // TPC \r
+ fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE); // to get chi2 and ncls of kTPCin\r
+ fEsdTrackCuts->SetMinNClustersTPC(70);\r
+ //fEsdTrackCuts->SetMinNClustersTPC(90); // ***** TMP *****\r
+ fEsdTrackCuts->SetMaxChi2PerClusterTPC(4);\r
+ fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);\r
+ fEsdTrackCuts->SetRequireTPCRefit(kTRUE);\r
+ // ITS\r
+ fEsdTrackCuts->SetRequireITSRefit(kTRUE);\r
+ fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);\r
+ fEsdTrackCuts->SetMaxDCAToVertexZ(2); // new for pile-up !!!\r
+\r
+ // Add SPD requirement \r
+ fEsdTrackCutsSPD = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");\r
+ fEsdTrackCutsSPD->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);\r
+ \r
+ // Add SDD requirement \r
+ fEsdTrackCutsSDD = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");\r
+ fEsdTrackCutsSDD->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);\r
+ \r
+ // Add DCA cuts \r
+ fEsdTrackCutsDCA = new AliESDtrackCuts("DCA", "pT dependent DCA cut");\r
+ // 7*(0.0050+0.0060/pt^0.9)\r
+ fEsdTrackCutsDCA->SetMaxDCAToVertexXYPtDep("0.0350+0.0420/pt^0.9");\r
+ \r
+ fEsdTrackCutsDCA->SetMaxDCAToVertexZ(1.e6);\r
+ fEsdTrackCutsDCA->SetDCAToVertex2D(kFALSE);\r
+\r
+ // emulates filterbit when getting leading-track from ESD\r
+ fAnalyseUE->DefineESDCuts(16); // any number = standard ITS+TPC + (SPD or SDD)\r
}
//____________________________________________________________________
if (fDebug > 1) AliInfo("CreateOutPutData()");
- // Create pointer to AliAnalysisUE, a class implementing the main analysis algorithms
- fAnaUE = new AliAnalyseUE();
- if (!fAnaUE){
+ // Create pointer to AliAnalyseLeadingTrackUE, a class implementing the main analysis algorithms\r
+ fAnalyseUE = new AliAnalyseLeadingTrackUE();\r
+ if (!fAnalyseUE){\r
AliError("UE analysis class not initialized!!!");
return;
}
- // Create pointer to AliHistogramsUE, a class handling histogram creation/filling
- fHistosUE = new AliHistogramsUE();
- if (!fHistosUE){
- AliError("UE histograms class not initialized!!!");
- return;
- }
-
+ \r
// Create list of output histograms
if (fListOfHistos != NULL){
delete fListOfHistos;
fListOfHistos->SetOwner(kTRUE);
}
+ // Create CF container\r
+ CreateContainer();\r
+ // number of events\r
+ fhEntries = new TH1F("fhEntries","Entries",1,0.,2.); \r
+ fListOfHistos->Add(fhEntries);\r
+ // tracks contributing to the vertex\r
+ fhVertexContributors = new TH1F("fhVertexContributors","Tracks in vertex",51, -0.5, 50.5);\r
+ fhVertexContributors->GetXaxis()->SetTitle("# tracks in vertex");\r
+ fhVertexContributors->GetYaxis()->SetTitle("Entries");\r
+ fListOfHistos->Add(fhVertexContributors);\r
+ // vertex resolution\r
+ fhVertexReso = new TH3F("fhVertexReso","Vertex resolution",51,-0.5,50.5,100,0.,0.05,100.,0.,0.1); \r
+ fhVertexContributors->GetXaxis()->SetTitle("# tracks in vertex");\r
+ fhVertexContributors->GetYaxis()->SetTitle("Resolution XY (cm)");\r
+ fhVertexContributors->GetZaxis()->SetTitle("Resolution Z (cm)");\r
+ fListOfHistos->Add(fhVertexReso);\r
- //Initialize output histograms
- fHistosUE->CreateHistogramsCorrections(fListOfHistos,fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, fTrackEtaCut);
+ // number of fake tracks\r
+ fhFakes = new TH1F("fhFakes","Fraction of fake tracks",5,-0.5,4.5);\r
+ fhFakes->GetXaxis()->SetBinLabel(1,"No MC");\r
+ fhFakes->GetXaxis()->SetBinLabel(2,"Unique MC primary");\r
+ fhFakes->GetXaxis()->SetBinLabel(3,"Unique MC secondary");\r
+ fhFakes->GetXaxis()->SetBinLabel(4,"Multiple MC");\r
+ fListOfHistos->Add(fhFakes);\r
+ \r
+ //pT distributions\r
+ fhPtMCAll = new TH1F("fhPtMCAll","All MC particles reconstructed",100,0., 20.);\r
+ fListOfHistos->Add(fhPtMCAll);\r
+ fhPtMCPrim = new TH1F("fhPtMCPrim","Primary MC particles reconstructed",100,0., 20.);\r
+ fListOfHistos->Add(fhPtMCPrim);\r
+ fhPtMCSec = new TH1F("fhPtMCSec","Secondary MC particles reconstructed",100,0., 20.);\r
+ fListOfHistos->Add(fhPtMCSec);\r
+ fhPtMCPrimFake = new TH1F("fhPtMCPrimFake","Fake primary MC particles reconstructed",100,0., 20.);\r
+ fListOfHistos->Add(fhPtMCPrimFake);\r
+ fhPtMCSecFake = new TH1F("fhPtMCSecFake","Fake secondary MC particles reconstructed",100,0., 20.);\r
+ fListOfHistos->Add(fhPtMCSecFake);\r
+\r
+\r
+ // Add task configuration to output list \r
AddSettingsTree();
-
- //Configure the CF manager
- if (fCFManager != NULL){
- delete fCFManager;
- fCFManager = NULL;
- }
- if (!fCFManager){
- fCFManager = new AliCFManager();
- }
- fHistosUE->CreateCorrectionsContainer(fCFManager,fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, fTrackEtaCut,fJet1EtaCut);
+ \r
//Post outputs
PostData(0,fListOfHistos);
- PostData(1,fCFManager->GetEventContainer());
-
+ \r
}
//____________________________________________________________________
void AliAnalysisTaskCorrectionsUE::Exec(Option_t */*option*/)
{
-
- Bool_t flag=kTRUE;
-
- //Determine the corrections
+ \r
+ // Get MC event\r
if (fMcHandler){
fMcEvent = fMcHandler->MCEvent();
if ( fDebug > 3 ) AliInfo( " Processing MC event..." );
- }else{
- if ( fDebug > 3 ) AliInfo( " Processing DATA event..." );
+ if (fMode && !fMcEvent) return;\r
}
-
- flag = EvaluateCorrections();
-
+ // Do the analysis\r
+ AnalyseCorrectionMode();\r
+
+ \r
+ PostData(0,fListOfHistos);\r
+
+}\r
+\r
+//____________________________________________________________________\r
+void AliAnalysisTaskCorrectionsUE::Terminate(Option_t */*option*/)\r
+{\r
+ // Terminate analysis\r
- if (flag){ //executed if event passes trigger+vertex selection
- //Fetch the pythia header info and get the trials
- Float_t nTrials = 1;
- if (fMcHandler && fMcEvent) {
- AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(fMcEvent);
- if(pythiaGenHeader) nTrials = pythiaGenHeader->Trials();
+ if (fDebug >1) AliAnalysisHelperJetTasks::PrintDirectorySize("PWG4_JetTasksOutput.root");\r
+ \r
+ if (!gROOT->IsBatch()){\r
+ fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));\r
+ if (!fListOfHistos){\r
+ AliError("Histogram List is not available");\r
+ return;\r
}
- fHistosUE->GetTrials()->Fill("#sum{ntrials}",fAvgTrials);
- PostData(0,fListOfHistos);
- PostData(1,fCFManager->GetEventContainer());
- }
+\r
+ } else {\r
+ AliInfo(" Batch mode, not histograms will be shown...");\r
+ }\r
+\r
+ \r
}
+\r
+/******************** ANALYSIS METHODS *****************************/\r
+\r
//____________________________________________________________________
void AliAnalysisTaskCorrectionsUE::AddSettingsTree()
{
//Write settings to output list
- TTree *settingsTree = new TTree("UECorrectionsSettings","Analysis Settings in UE corrections");
- settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
- settingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
- settingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D");
- settingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D");
- settingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D");
- settingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D");
+ TTree *settingsTree = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");\r
+ settingsTree->Branch("fnTracksVertex", &fnTracksVertex, "TracksInVertex/D");\r
+ settingsTree->Branch("fZVertex", &fZVertex, "VertexZCut/D");\r
settingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D");
settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
- settingsTree->Branch("fAnaType", &fAnaType, "Ana/I");
- settingsTree->Branch("fRegionType", &fRegionType,"Reg/I");
- settingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I");
- settingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O");
- settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
- settingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O");
- settingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O");
settingsTree->Fill();
fListOfHistos->Add(settingsTree);
-}
-
+} \r
//____________________________________________________________________
-Bool_t AliAnalysisTaskCorrectionsUE::EvaluateCorrections()
+void AliAnalysisTaskCorrectionsUE::AnalyseCorrectionMode()\r
{
-
- /////////////////////////////////////////////////////////////////////////
- //
- // EVENT SELECTION:
- // CF containers are filled to get the number of entries after every cut.
- // *** Cuts: ***
- // 0 - triggered
- // 1 - physics selection
- // 2 - vertex selection
- // 3 - event topology
- // 4 - leading track pT cut
- // 5 - leading track correctly identified
- // *** Variables: ***
- // 0 - leading track pT (reco)
- // 1 - leading track eta (reco)
- // 2 - process type:
- // 1: non diffractive
- // 2: double diffractive
- // 4: single diffractive
- // 3 - leading track pT (true)
- // 4 - leading track eta (true)
- // 5 - delta eta (reco-true)
- // 6 - delta phi (reco-true)
- // 7 - radius (reco-true)
- //
- // TRACK-LEVEL CORRECTIONS:
- // Fill histograms similar to AliAnalysisTaskUE.
- //
- /////////////////////////////////////////////////////////////////////////
- Double_t containerInput[8];// relevant variables (see above)
-
- //PROCESS TYPE (ND,SD,DD)
- AliAnalysisHelperJetTasks::MCProcessType eventId = AliAnalysisHelperJetTasks::kInvalidProcess;
- if (fMcHandler && fMcEvent) {
- AliGenEventHeader* genHeader = fMcEvent->GenEventHeader();
- eventId = AliAnalysisHelperJetTasks::GetPythiaEventProcessType(genHeader,kFALSE);
- if (eventId<0){
- eventId = AliAnalysisHelperJetTasks::GetDPMjetEventProcessType(genHeader,kFALSE);
+ // Analyse the event\r
+ Int_t labelMC = -1;\r
+ if (fMcHandler && fMcEvent){\r
+ // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex\r
+ if (!fAnalyseUE->VertexSelection(fMcEvent, 0, fZVertex)) \r
+ return; \r
+ // Get MC-true leading particle \r
+ TObjArray *ltMC = (TObjArray*)fAnalyseUE->FindLeadingObjects(fMcEvent);\r
+ AliVParticle* leadingMC = 0;\r
+ if (ltMC){\r
+ leadingMC = (AliVParticle*) ltMC->At(0);\r
}
- if (eventId<0 && fDebug>1)AliInfo("No Pythia or Phojet Header retrived!");
- }else if (fDebug>1) AliInfo("No MC handler or Event!");
-
- //Initialize container inputs
- for (Int_t i =0; i<8; i++){
- containerInput[i]=-999.;
+ if (!leadingMC)return; \r
+ labelMC = TMath::Abs(leadingMC->GetLabel());\r
+ }\r
+
+ // Trigger selection ************************************************\r
+ if (!fAnalyseUE->TriggerSelection(fESDHandler)) return;\r
+ \r
+ // PILEUP-CUT ****** NEW !!!!!!!!! **************\r
+ Bool_t select = kTRUE;\r
+ //select = AliAnalysisHelperJetTasks::TestSelectInfo(AliAnalysisHelperJetTasks::kIsPileUp);\r
+ if (! select) return;\r
+ \r
+ // Vertex selection *************************************************\r
+ \r
+ if(!fAnalyseUE->VertexSelection(fESDEvent, 0, fZVertex)) return;\r
+ AliESDVertex* vertex = (AliESDVertex*)fESDEvent->GetPrimaryVertex();\r
+ Int_t nvtx = vertex->GetNContributors();\r
+ fhVertexContributors->Fill(nvtx);\r
+ if (fMcHandler){\r
+ AliVVertex *vertexMC = (AliVVertex*)fMcEvent->GetPrimaryVertex();\r
+ if (vertexMC){\r
+ Double_t diffx = vertexMC->GetX()-vertex->GetX();\r
+ Double_t diffy = vertexMC->GetY()-vertex->GetY();\r
+ Double_t diffxy = TMath::Sqrt(TMath::Power(diffx,2)+TMath::Power(diffy,2));\r
+ //Double_t diffxy = TMath::Abs(diffx); // **** TMP ****\r
+ //Double_t diffxy = diffy;\r
+ Double_t diffz = TMath::Abs(vertexMC->GetZ()-vertex->GetZ());\r
+ \r
+ fhVertexReso->Fill(nvtx,diffxy,diffz);\r
+ }else if (fDebug>1)Printf("******* NO MC VERTEX ********");\r
}
+\r
+ if(!fAnalyseUE->VertexSelection(fESDEvent, fnTracksVertex, fZVertex)) return;\r
+\r
+ // Get Reconstructed leading particle *******************************\r
+ TObjArray *ltRECO = fAnalyseUE->FindLeadingObjects(fESDEvent);\r
+ if (!ltRECO){\r
+ delete[] ltRECO;\r
+ return;\r
+ }\r
+ Int_t labelReco= TMath::Abs(((AliVParticle*)ltRECO->At(0))->GetLabel());\r
- //Assign process type
- if (eventId == 1 ) containerInput[2]=1.; //Non diffractive
- if (eventId == 2 ) containerInput[2]=2.; //Double diffractive
- if (eventId == 4 ) containerInput[2]=4.; //Single diffractive
- // Execute analysis for current event ******************************
+ // Loop on tracks\r
+ Int_t nTracks = fESDEvent->GetNumberOfTracks();\r
+ if (!nTracks)return;\r
+ // count accepted events\r
+ fhEntries->Fill(1.);\r
+\r
+ Int_t npart=0;\r
+ Bool_t *labelsArray = 0; \r
+ if (fMcHandler){\r
+ npart = fMcEvent->GetNumberOfTracks();\r
+ labelsArray = new Bool_t[npart];\r
+ for(Int_t j = 0; j<npart; j++){\r
+ labelsArray[j] = kFALSE;\r
+ }\r
+ }\r
+\r
+ for (Int_t i = 0; i < nTracks; i++){\r
+ AliESDtrack *track = fESDEvent->GetTrack(i);\r
+ // only charged\r
+ if (!track || !track->Charge())continue;\r
+ // apply cuts\r
+ Bool_t cut = fEsdTrackCuts->AcceptTrack(track);\r
+ Bool_t cutSPD = fEsdTrackCutsSPD->AcceptTrack(track);\r
+ Bool_t cutSDD = fEsdTrackCutsSDD->AcceptTrack(track);\r
+ Bool_t cutDCA = fEsdTrackCutsDCA->AcceptTrack(track);\r
+ \r
+ //Exclude the MC leading track\r
+ Double_t matchLeading = 1.;//no match\r
+ if (fMcHandler){\r
+ if (TMath::Abs(track->GetLabel())==labelMC) matchLeading=0.; //match MC leading\r
+ if (TMath::Abs(track->GetLabel())==labelReco) {\r
+ matchLeading=2.;//match RECO leading\r
+ if (labelMC == labelReco)matchLeading = 3.; // match both (mc = reco leading)\r
+ }\r
+ }\r
+ // Fill step0 (all tracks) \r
+ FillContainer(track, 0,kFALSE,matchLeading); \r
+ // Fill step1 (no SPD cluster requirement - no DCA cut )\r
+ if ( cut ) FillContainer(track,1,kFALSE,matchLeading);\r
+ // Fill step2\r
+ if ( cut && cutDCA && (cutSPD || cutSDD)) FillContainer(track,2,kFALSE,matchLeading);\r
+ // Fill step3-step4-step5 \r
+ if ( cut && cutDCA && !cutSPD && !cutSDD) FillContainer(track,3,kTRUE,matchLeading);\r
+ if ( cut && cutDCA && cutSPD) FillContainer(track,4,kTRUE,matchLeading);\r
+ if ( cut && cutDCA && !cutSPD && cutSDD) FillContainer(track,5,kTRUE,matchLeading);\r
+ // Fill step 6 - temporary just to define the standard track cuts \r
+ if ( cut && cutSPD ) FillContainer(track,6,kFALSE,matchLeading);\r
+ //if ( cut && cutSPD ) FillContainer(track,6,kTRUE,matchLeading); // ***** TMP *****\r
+ if ( cut && (cutSPD || cutSDD) ) FillContainer(track,7,kFALSE,matchLeading);\r
+ // Study contamination form fakes\r
+ if (fMcHandler){\r
+ \r
+ // consider standard ITS+TPC cuts without SPD cluster requirement\r
+ if (cut && cutDCA){\r
+ //check if it points back to any MC\r
+ Int_t label = TMath::Abs(track->GetLabel());\r
+ AliVParticle *part = (AliVParticle*)fMcEvent->GetTrack(label); \r
+ if (!part){\r
+ fhFakes->Fill(0.);\r
+ //Printf("*************** NO MC PARTICLE ************************");\r
+ continue;\r
+ }\r
+ \r
+ fhPtMCAll->Fill(part->Pt()); \r
+ // Check if label is not already in array\r
+ if (!labelsArray[label]){\r
+ labelsArray[label]= kTRUE;\r
+ if (fMcEvent->IsPhysicalPrimary(label)){\r
+ fhFakes->Fill(1.);\r
+ fhPtMCPrim->Fill(part->Pt());\r
+ }else{\r
+ fhFakes->Fill(2.);\r
+ fhPtMCSec->Fill(part->Pt());\r
+ }\r
+ }else{\r
+ fhFakes->Fill(3.);\r
+ if (fMcEvent->IsPhysicalPrimary(label))fhPtMCPrimFake->Fill(part->Pt());\r
+ else fhPtMCSecFake->Fill(part->Pt());\r
+ }\r
+ } \r
+ }\r
+ } // end loop on tracks\r
+ if(labelsArray)\r
+ delete[] labelsArray;\r
+\r
+}\r
+\r
+\r
+//____________________________________________________________________\r
+void AliAnalysisTaskCorrectionsUE::CreateContainer()\r
+{\r
- // Get jets and order by pT
- TVector3 jetVect[3];
- *jetVect = fAnaUE->GetOrderedClusters(fAODBranch, fUseChPartJet, fChJetPtMin );
+ // Create the output CF container\r
+ // relevant variables \r
+ UInt_t iprim = 0; // 0: primaries, 1: secondaries from strangness, 2: secondaries form material \r
+ UInt_t iptmc = 1;\r
+ UInt_t ipt = 2;\r
+ UInt_t ieta = 3;\r
+ UInt_t idcaxy = 4;\r
+ UInt_t idcaz = 5;\r
+ UInt_t imatch = 6;\r
+ UInt_t icharge = 7;\r
+ // set-up the grid\r
+ UInt_t nstep = 8;\r
+ const Int_t nvar = 8;\r
+ const Int_t nbin0 = 5; // prim\r
+ const Int_t nbin1 = 20; // pt resolution\r
+ const Int_t nbin2 = 39; // pt \r
+ const Int_t nbin3 = 20; // eta\r
+ const Int_t nbin4 = 100; // dca xy\r
+ const Int_t nbin5 = 100; // dca z\r
+ const Int_t nbin6 = 4; // matching with leading track\r
+ const Int_t nbin7 = 2;\r
+\r
+ // array for the number of bins in each dimension\r
+ Int_t iBin[nvar];\r
+ iBin[0] = nbin0;\r
+ iBin[1] = nbin1;\r
+ iBin[2] = nbin2;\r
+ iBin[3] = nbin3;\r
+ iBin[4] = nbin4;\r
+ iBin[5] = nbin5;\r
+ iBin[6] = nbin6;\r
+ iBin[7] = nbin7;\r
- //now define leading track pT and eta
- containerInput[0]=jetVect[0].Pt();
- containerInput[1]=jetVect[0].Eta();
-
- fCFManager->GetEventContainer()->Fill(containerInput,kCFStepTriggered); //fill CF container
+ // primaries\r
+ Double_t primBins[7] = {-0.5,0.5,1.5,2.5,3.5,4.5,5.5};\r
+ // matching with leading-track\r
+ Double_t matchBins[5] = {-0.5,0.5,1.5,2.5,3.5};\r
+
+ // pT resolution\r
+ Double_t resoBins[nbin1+1];\r
+ for (Int_t i=0; i<=nbin1; i++)\r
+ resoBins[i] = -1.0 + 0.1 * i;\r
- // Physics selection ************************************************
- fInputHandler = (AliInputEventHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
- if (!fAnaUE->TriggerSelection(fInputHandler)) return kFALSE;
+ // pT\r
+ Double_t ptBins[nbin2+1] = {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};\r
- fCFManager->GetEventContainer()->Fill(containerInput, kCFStepPhysSelect); //fill CF container
+ // eta\r
+ Double_t etaBins[nbin3+1];\r
+ for (Int_t i=0; i<=nbin3; i++)\r
+ etaBins[i] = -1.0 + 0.1 * i;\r
- // Event selection (vertex) *****************************************
+ // dca xy\r
+ Double_t dcaxyBins[nbin4+1];\r
+ for (Int_t i=0; i<=nbin4; i++)\r
+ dcaxyBins[i] = -1.+0.02 * i;\r
- if(!fAnaUE->VertexSelection(fAOD,fnTracksVertex,fZVertex)) return kFALSE;
- //if(!fAnaUE->VertexSelectionOld(fAOD)) return kFALSE; // temporary to compare with old task and to have same cuts for MC !!!
-
- fCFManager->GetEventContainer()->Fill(containerInput, kCFStepVertexSelect); //fill CF container
-
- // Select events according to analysis type *************************
- // (in the leading track analysis it should not happen that there are no "jets")
- if( jetVect[0].Pt() < 0. ) {
- if( fDebug > 1 ) AliInfo("\n Skipping Event, not jet found...");
- return kTRUE;
- } else {
- if (fDebug >1 ) AliInfo(Form("\n Pt Leading Jet = %6.1f eta=%5.3f ", jetVect[0].Pt(), jetVect[0].Eta() ));
- }
+ // dca z\r
+ Double_t dcazBins[nbin5+1];\r
+ for (Int_t i=0; i<=nbin5; i++)\r
+ dcazBins[i] = -5.0 + 0.1 * i;\r
+ \r
+ // charge \r
+ Double_t chargeBins[nbin7+1] = {-1.,0.,1.};\r
+
+ // create container\r
+ // set variables\r
+ fOutCFcont = new AliCFContainer("fOutCFcont","Output Container",nstep,nvar,iBin);\r
+ fOutCFcont->SetBinLimits(iprim,primBins);\r
+ fOutCFcont->SetVarTitle(iprim, "Particle type");\r
+ fOutCFcont->SetBinLimits(iptmc,resoBins);\r
+ fOutCFcont->SetVarTitle(iptmc, "#Delta p_{T} (DATA-MC) (GeV/c)");\r
+ fOutCFcont->SetBinLimits(ipt,ptBins);\r
+ fOutCFcont->SetVarTitle(ipt, "p_{T} (GeV/c)");\r
+ fOutCFcont->SetBinLimits(ieta,etaBins);\r
+ fOutCFcont->SetVarTitle(ieta, "#eta");\r
+ fOutCFcont->SetBinLimits(idcaxy,dcaxyBins);\r
+ fOutCFcont->SetVarTitle(idcaxy, " DCA_{XY} (cm)");\r
+ fOutCFcont->SetBinLimits(idcaz,dcazBins);\r
+ fOutCFcont->SetVarTitle(idcaz, " DCA_{Z} (cm)");\r
+ fOutCFcont->SetBinLimits(imatch,matchBins);\r
+ fOutCFcont->SetVarTitle(imatch, "Matching with leading-track");\r
+ fOutCFcont->SetBinLimits(icharge,chargeBins);\r
+ fOutCFcont->SetVarTitle(icharge, "Charge");\r
+\r
+ // set steps\r
+ fOutCFcont->SetStepTitle(0,"all tracks");\r
+ fOutCFcont->SetStepTitle(1,"ITS+TPC cuts (no SPD cluster requirement and DCA cut)");\r
+ fOutCFcont->SetStepTitle(2,"add DCA cut");\r
+ fOutCFcont->SetStepTitle(3,"NO SPD cluster, NO SDD cluster in first layer");\r
+ fOutCFcont->SetStepTitle(4,"YES SPD cluster, NO SDD cluster in first layer");\r
+ fOutCFcont->SetStepTitle(5,"NO SPD cluster, YES SDD cluster in first layer");\r
+ fOutCFcont->SetStepTitle(6,"ITS+TPC cuts - no DCA cut - SPD cut");\r
+ fOutCFcont->SetStepTitle(7,"ITS+TPC cuts - no DCA cut - SPD or SDD cut");\r
+ fListOfHistos->Add(fOutCFcont);\r
+\r
+}\r
+\r
+\r
+void AliAnalysisTaskCorrectionsUE::FillContainer(AliESDtrack *track, Int_t step,Bool_t mcvertex, Double_t matchLeading)\r
+{\r
+\r
+ // Fill the CF container\r
+\r
+ Double_t vars[8];\r
+ Double_t prim = -1.;\r
+ if (track->Charge() > 0.) vars[7] = 0.5;\r
+ else vars[7] = -0.5;\r
- if ( ! (fAnaUE->AnaTypeSelection(jetVect))) return kTRUE;
-
-
- // For the events selected check the real MC leading particle
- // (only when running on MC)
- TVector3 jetVectMC[3];
- if (fMcEvent){
- *jetVectMC = fAnaUE->GetLeadingTracksMC(fMcEvent);
- if (fAnaUE->AnaTypeSelection(jetVectMC) ){
- //now define leading track pT and eta
- containerInput[3]=jetVectMC[0].Pt();
- containerInput[4]=jetVectMC[0].Eta();
- //Check distance between real and reco leading-particle
- containerInput[5] = jetVect[0].Eta()-jetVectMC[0].Eta();
- containerInput[6] = TMath::Abs(jetVect[0].Phi()-jetVectMC[0].Phi());
- if (containerInput[6] >= 2.*TMath::Pi())containerInput[6] -= 2.*TMath::Pi();
- containerInput[7] = sqrt(TMath::Power(containerInput[5],2.) + TMath::Power(containerInput[6],2.));
- }
- }
-
- fCFManager->GetEventContainer()->Fill(containerInput, kCFStepAnaTopology); //fill CF container
-
- //Cut on the leading track pT
- if (!(jetVect[0].Pt() >= 1.)) return kTRUE;
- // (only when running on MC)
- if (fMcEvent){
- if (!(jetVectMC[0].Pt()>= 1.)){
- containerInput[3]=containerInput[4]=-999.;
- containerInput[5]=containerInput[6]=containerInput[7]=-999.;
+ if (fMcHandler){\r
+ // determine if points back to a primary\r
+ Int_t label = TMath::Abs(track->GetLabel());\r
+ \r
+ AliMCParticle *part = (AliMCParticle*)fMcEvent->GetTrack(label); \r
+ for (Int_t i=0; i<=6;i++) vars[i] = -999.;\r
+ if (part) { //PRIMARY\r
+ if (fMcEvent->IsPhysicalPrimary(label)){\r
+ prim = 0.;\r
+ }else { //SECONDARY\r
+ // decide if strange\r
+ Int_t labelm = TMath::Abs(part->GetMother());\r
+ AliMCParticle *mother = (AliMCParticle*)fMcEvent->GetTrack(labelm); \r
+ Int_t code = mother->PdgCode();\r
+ Int_t mfl = Int_t (code/ TMath::Power(10, Int_t(TMath::Log10(code))));\r
+ if (mfl == 3) prim = 1.;\r
+ else{ \r
+ //Printf("***** PROCESS : %d",part->Particle()->GetUniqueID()); \r
+ if (TMath::Abs(code) == 211) prim = 2.; // charged pion decay\r
+ else if (part->Particle()->GetUniqueID() == 13 )prim = 3.; // hadronic interactions\r
+ else if (part->Particle()->GetUniqueID() == 5 )prim = 4.; // photon conversions\r
+ else prim = 5.; // other?\r
+ }\r
+ }\r
+ vars[1]= part->Pt()-track->Pt();\r
+ // In step 2 fill MC pT for contamination study\r
+ if (step == 2)vars[2]=part->Pt();\r
}
- }
-
- fCFManager->GetEventContainer()->Fill(containerInput, kCFStepLtPtCut1); //fill CF container
-
- // Check if leading track is correctly identified
- // (only when running on MC)
- if (fMcEvent){
- Int_t labelLt = fAnaUE->GetLtLabel();
- Int_t labelLtMC = fAnaUE->GetLtMCLabel();
- if (labelLt == labelLtMC){
- fCFManager->GetEventContainer()->Fill(containerInput, kCFStepLtCorrect); //fill CF container
- }
- }
-
- // For track efficiency and contamination
- // (only when running on MC)
- if (fMcEvent){
- //Run once on reco...
- fAnaUE->FindMaxMinRegions( jetVect, fConePosition, 0, 0 );//normal track cut
- fAnaUE->FindMaxMinRegions( jetVect, fConePosition, 0, 1 );//for efficiency: cut on pT and eta are set on MC true
- //and once on MC true
- fAnaUE->FindMaxMinRegions( jetVect, fConePosition, 1, 0 );
- }else{
- // For d0 distribution, runs only on real data
- fAnaUE->FindMaxMinRegions( jetVect, fConePosition, 0, 0 );//run on real data
- }
-
-
- return kTRUE;
+ }\r
+ vars[0]=prim;\r
+
+ if (step != 2)vars[2]=track->Pt();\r
+ vars[3]=track->Eta();\r
+
+ Bool_t dcaControlFlag = kFALSE;\r
+ if (mcvertex && fMcHandler){\r
+ // we want DCA w.r.t. MC vertex\r
+ AliVVertex *vtxMC = (AliVVertex*)fMcEvent->GetPrimaryVertex();\r
+ dcaControlFlag = track->RelateToVertex((AliESDVertex*)vtxMC,(Double_t)fESDEvent->GetMagneticField(),10000.);\r
+ }else{\r
+ AliESDVertex* vertex = (AliESDVertex*)fESDEvent->GetPrimaryVertex();\r
+ dcaControlFlag = track->RelateToVertex(vertex,(Double_t)fESDEvent->GetMagneticField(),10000.);\r
+ }\r
+ if (dcaControlFlag){\r
+ Float_t dca[2];\r
+ Float_t dcaCov[2];\r
+ track->GetImpactParameters(dca,dcaCov);\r
+ vars[4]=dca[0];\r
+ vars[5]=dca[1];\r
+ } \r
+
+
+ vars[6]= matchLeading;\r
+ fOutCFcont->Fill(vars,step);\r
+\r
}
+\r
//____________________________________________________________________
AliAnalysisTaskCorrectionsUE* AliAnalysisTaskCorrectionsUE::Instance()
{
}
}
-//____________________________________________________________________
-void AliAnalysisTaskCorrectionsUE::Terminate(Option_t */*option*/)
-{
- // Terminate analysis
-
- if (fDebug >1) AliAnalysisHelperJetTasks::PrintDirectorySize("PWG4_JetTasksOutput.root");
-
- if (!gROOT->IsBatch()){
- fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
- if (!fListOfHistos){
- AliError("Histogram List is not available");
- return;
- }
-
- } else {
- AliInfo(" Batch mode, not histograms will be shown...");
- }
-
-
-}
-
-//____________________________________________________________________
-void AliAnalysisTaskCorrectionsUE::WriteSettings()
-{
-
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff