From 8838ab7ab4663d7ad4117ce1b692143175fa0d4e Mon Sep 17 00:00:00 2001 From: morsch Date: Tue, 12 May 2009 09:32:11 +0000 Subject: [PATCH] Updates needed for full jet reconstruction (charged + emcal) [Magali Estienne] --- JETAN/AliAnalysisTaskJets.cxx | 157 +- JETAN/AliAnalysisTaskJets.h | 16 +- JETAN/AliFastJetFinder.cxx | 282 +++- JETAN/AliFastJetFinder.h | 13 +- JETAN/AliFastJetHeaderV1.cxx | 89 + JETAN/AliFastJetHeaderV1.h | 97 ++ JETAN/AliJet.cxx | 67 +- JETAN/AliJet.h | 78 +- JETAN/AliJetESDReader.cxx | 174 +- JETAN/AliJetESDReader.h | 62 +- JETAN/AliJetFillUnitArrayEMCalDigits.cxx | 438 +++-- JETAN/AliJetFillUnitArrayEMCalDigits.h | 67 +- JETAN/AliJetFillUnitArrayTracks.cxx | 476 ++++-- JETAN/AliJetFillUnitArrayTracks.h | 112 +- JETAN/AliJetFinder.cxx | 130 +- JETAN/AliJetFinder.h | 57 +- JETAN/AliJetGrid.cxx | 226 +-- JETAN/AliJetGrid.h | 25 +- JETAN/AliJetReader.cxx | 57 +- JETAN/AliJetReader.h | 74 +- JETAN/AliJetReaderHeader.cxx | 2 + JETAN/AliJetReaderHeader.h | 9 +- JETAN/AliJetUnitArray.cxx | 200 ++- JETAN/AliJetUnitArray.h | 109 +- JETAN/AliSISConeJetFinder.cxx | 294 +++- JETAN/AliSISConeJetFinder.h | 9 +- JETAN/AliSISConeJetHeader.cxx | 43 +- JETAN/AliSISConeJetHeader.h | 36 +- JETAN/AliUA1JetFinderV2.cxx | 1939 ++++++++++++++++------ JETAN/AliUA1JetFinderV2.h | 29 +- 30 files changed, 3836 insertions(+), 1531 deletions(-) create mode 100644 JETAN/AliFastJetHeaderV1.cxx create mode 100644 JETAN/AliFastJetHeaderV1.h diff --git a/JETAN/AliAnalysisTaskJets.cxx b/JETAN/AliAnalysisTaskJets.cxx index 4e755a9d426..b34436e68f7 100644 --- a/JETAN/AliAnalysisTaskJets.cxx +++ b/JETAN/AliAnalysisTaskJets.cxx @@ -14,7 +14,7 @@ **************************************************************************/ /* $Id$ */ - +#include #include #include #include @@ -26,6 +26,8 @@ #include "AliAnalysisManager.h" #include "AliJetFinder.h" #include "AliJetHeader.h" +#include "AliJetReader.h" +#include "AliJetReaderHeader.h" #include "AliJetHistos.h" #include "AliESDEvent.h" #include "AliESD.h" @@ -42,114 +44,143 @@ ClassImp(AliAnalysisTaskJets) //////////////////////////////////////////////////////////////////////// AliAnalysisTaskJets::AliAnalysisTaskJets(): - AliAnalysisTaskSE(), - fConfigFile("ConfigJetAnalysis.C"), - fNonStdBranch(""), + AliAnalysisTaskSE(), + fConfigFile("ConfigJetAnalysis.C"), + fNonStdBranch(""), fJetFinder(0x0), fHistos(0x0), - fListOfHistos(0x0) + fListOfHistos(0x0), + fChain(0x0), + fOpt(0) { // Default constructor } AliAnalysisTaskJets::AliAnalysisTaskJets(const char* name): - AliAnalysisTaskSE(name), - fConfigFile("ConfigJetAnalysis.C"), - fNonStdBranch(""), - fJetFinder(0x0), - fHistos(0x0), - fListOfHistos(0x0) + AliAnalysisTaskSE(name), + fConfigFile("ConfigJetAnalysis.C"), + fNonStdBranch(""), + fJetFinder(0x0), + fHistos(0x0), + fListOfHistos(0x0), + fChain(0x0), + fOpt(0) { // Default constructor - DefineOutput(1, TList::Class()); + DefineOutput(1, TList::Class()); } -void AliAnalysisTaskJets::UserCreateOutputObjects() +AliAnalysisTaskJets::AliAnalysisTaskJets(const char* name, TChain* chain): + AliAnalysisTaskSE(name), + fConfigFile("ConfigJetAnalysis.C"), + fNonStdBranch(""), + fJetFinder(0x0), + fHistos(0x0), + fListOfHistos(0x0), + fChain(chain), + fOpt(0) { -// Create the output container -// - if (fDebug > 1) printf("AnalysisTaskJets::CreateOutPutData() \n"); - - + // Default constructor + DefineOutput(1, TList::Class()); +} - if(fNonStdBranch.Length()==0){ +//---------------------------------------------------------------- +void AliAnalysisTaskJets::UserCreateOutputObjects() +{ + // Create the output container + // + if (fDebug > 1) printf("AnalysisTaskJets::CreateOutPutData() \n"); + + if(fNonStdBranch.Length()==0) + { // Connec default AOD to jet finder fJetFinder->ConnectAOD(AODEvent()); } - else{ + else + { // Create a new branch for jets... - // how is this is reset cleared in the UserExec.... + // how is this reset? -> cleared in the UserExec.... // Can this be handled by the framework? TClonesArray *tca = new TClonesArray("AliAODJet", 0); tca->SetName(fNonStdBranch); AddAODBranch("TClonesArray",&tca); - fJetFinder->ConnectAODNonStd(AODEvent(),fNonStdBranch.Data()); + fJetFinder->ConnectAODNonStd(AODEvent(), fNonStdBranch.Data()); + } + + // Histograms + OpenFile(1); + fListOfHistos = new TList(); + fHistos = new AliJetHistos(); + fHistos->AddHistosToList(fListOfHistos); + + // Add the JetFinderInformation to the Outputlist + AliJetHeader *fH = fJetFinder->GetHeader(); + + // Compose a characteristic output name + // with the name of the output branch + if(fH) { + if(fNonStdBranch.Length()==0) { + fH->SetName("AliJetHeader_jets"); } - // Histograms - OpenFile(1); - fListOfHistos = new TList(); - fHistos = new AliJetHistos(); - fHistos->AddHistosToList(fListOfHistos); - - // Add the JetFinderInforamtion to the Outputlist - AliJetHeader *fH = fJetFinder->GetHeader(); - // Compose a characteristic output name - // with the name of the output branch - if(fH){ - if(fNonStdBranch.Length()==0){ - fH->SetName("AliJetHeader_jets"); - } - else{ - fH->SetName(Form("AliJetHeader_%s",fNonStdBranch.Data())); - } + else { + fH->SetName(Form("AliJetHeader_%s",fNonStdBranch.Data())); } - OutputTree()->GetUserInfo()->Add(fH); + } + OutputTree()->GetUserInfo()->Add(fH); } +//---------------------------------------------------------------- void AliAnalysisTaskJets::Init() { - // Initialization - if (fDebug > 1) printf("AnalysisTaskJets::Init() \n"); - - // Call configuration file - if (fConfigFile.Length()) { - gROOT->LoadMacro(fConfigFile); - fJetFinder = (AliJetFinder*) gInterpreter->ProcessLine("ConfigJetAnalysis()"); - } - // Initialise Jet Analysis - fJetFinder->Init(); - // Write header information to local file - fJetFinder->WriteHeaders(); + // Initialization + if (fDebug > 1) printf("AnalysisTaskJets::Init() \n"); + + // Call configuration file + if (fConfigFile.Length()) { + gROOT->LoadMacro(fConfigFile); + fJetFinder = (AliJetFinder*) gInterpreter->ProcessLine("ConfigJetAnalysis()"); + } + AliJetReaderHeader *header = (AliJetReaderHeader*)fJetFinder->GetReader()->GetReaderHeader(); + fOpt = header->GetDetector(); + + // Initialise Jet Analysis + if(fOpt == 0) fJetFinder->Init(); + else fJetFinder->InitTask(fChain); // V2 + + // Write header information to local file + fJetFinder->WriteHeaders(); } - - - - +//---------------------------------------------------------------- void AliAnalysisTaskJets::UserExec(Option_t */*option*/) { // Execute analysis for current event // - - // Fill control histos TClonesArray* jarray = 0; - if(fNonStdBranch.Length()==0){ + + if(fNonStdBranch.Length()==0) { jarray = AODEvent()->GetJets(); } - else{ - jarray = dynamic_cast(AODEvent()->FindListObject(fNonStdBranch.Data())); + else { + jarray = dynamic_cast(AODEvent()->FindListObject(fNonStdBranch.Data())); jarray->Delete(); // this is our responsibility, clear before filling again } fJetFinder->GetReader()->SetInputEvent(InputEvent(), AODEvent(), MCEvent()); - fJetFinder->ProcessEvent(); - + if(fOpt==0) fJetFinder->ProcessEvent(); + else fJetFinder->ProcessEvent2(); // V2 + + // Fill control histos fHistos->FillHistos(jarray); + // Post the data PostData(1, fListOfHistos); } + +//************************************************************* + void AliAnalysisTaskJets::Terminate(Option_t */*option*/) { // Terminate analysis diff --git a/JETAN/AliAnalysisTaskJets.h b/JETAN/AliAnalysisTaskJets.h index 6d94463f634..02135d72673 100644 --- a/JETAN/AliAnalysisTaskJets.h +++ b/JETAN/AliAnalysisTaskJets.h @@ -8,16 +8,17 @@ class AliJetFinder; class AliESDEvent; class TTree; +class TChain; class AliAODEvent; class AliJetHistos; - class AliAnalysisTaskJets : public AliAnalysisTaskSE { public: AliAnalysisTaskJets(); AliAnalysisTaskJets(const char* name); + AliAnalysisTaskJets(const char* name, TChain* chain); virtual ~AliAnalysisTaskJets() {;} // Implementation of interface methods virtual void UserCreateOutputObjects(); @@ -25,19 +26,22 @@ class AliAnalysisTaskJets : public AliAnalysisTaskSE virtual void LocalInit() {Init();} virtual void UserExec(Option_t *option); virtual void SetConfigFile(const char *c){fConfigFile = c;} - void SetJetFinder(AliJetFinder *finder) {fJetFinder = finder;} virtual void SetNonStdBranch(const char *c){fNonStdBranch = c;} virtual void Terminate(Option_t *option); + private: AliAnalysisTaskJets(const AliAnalysisTaskJets &det); AliAnalysisTaskJets &operator=(const AliAnalysisTaskJets &det); private: TString fConfigFile; // the name of the ConfigFile - TString fNonStdBranch; // the name of the non-std branch name - AliJetFinder* fJetFinder; // Pointer to the jet finder - AliJetHistos* fHistos; // Histogram manager class - TList* fListOfHistos; // Output list of histograms + TString fNonStdBranch; // the name of the non-std branch name//commented by syssy + AliJetFinder* fJetFinder; // Pointer to the jet finder + AliJetHistos* fHistos; // Histogram manager class + TList* fListOfHistos; // Output list of histograms + TChain* fChain; // Chain + Int_t fOpt; // Detector configuration used + ClassDef(AliAnalysisTaskJets, 3); // Analysis task for standard jet analysis }; diff --git a/JETAN/AliFastJetFinder.cxx b/JETAN/AliFastJetFinder.cxx index 833bfa5815e..69d480144ec 100644 --- a/JETAN/AliFastJetFinder.cxx +++ b/JETAN/AliFastJetFinder.cxx @@ -16,11 +16,14 @@ //--------------------------------------------------------------------- // FastJet v2.3.4 finder algorithm interface +// Last modification: Neutral cell energy included in the jet reconstruction +// +// Authors: Rafael.Diaz.Valdes@cern.ch +// Magali.estienne@subatech.in2p3.fr (neutral part + bg subtraction option) // -// Author: Rafael.Diaz.Valdes@cern.ch -// //--------------------------------------------------------------------- + #include #include #include @@ -31,11 +34,11 @@ #include #include "AliFastJetFinder.h" -#include "AliFastJetHeader.h" +#include "AliFastJetHeaderV1.h" #include "AliJetReaderHeader.h" #include "AliJetReader.h" #include "AliJet.h" - +#include "AliJetUnitArray.h" #include "fastjet/PseudoJet.hh" #include "fastjet/ClusterSequenceArea.hh" @@ -56,6 +59,8 @@ using namespace std; ClassImp(AliFastJetFinder) + + //____________________________________________________________________________ AliFastJetFinder::AliFastJetFinder(): @@ -75,10 +80,11 @@ AliFastJetFinder::~AliFastJetFinder() void AliFastJetFinder::FindJets() { - Bool_t debug = kFALSE; - //pick up fastjet header - AliFastJetHeader *header = (AliFastJetHeader*)fHeader; + AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; + Bool_t debug = header->GetDebug(); // debug option + Bool_t bgMode = header->GetBGMode(); // choose to subtract BG or not + Int_t fOpt = fReader->GetReaderHeader()->GetDetector(); // check if we are reading AOD jets TRefArray *refs = 0; @@ -88,26 +94,62 @@ void AliFastJetFinder::FindJets() // RUN ALGORITHM // read input particles ----------------------------- vector input_particles; - TClonesArray *lvArray = fReader->GetMomentumArray(); - if(lvArray == 0) { cout << "Could not get the momentum array" << endl; return; } - Int_t nIn = lvArray->GetEntries(); - if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; } - //Int_t nJets = 0; // n jets in this event - fJets->SetNinput(nIn) ; // number of input objects - Float_t px,py,pz,en; - // load input vectors - for(Int_t i = 0; i < nIn; i++){ // loop for all input particles - TLorentzVector *lv = (TLorentzVector*) lvArray->At(i); - px = lv->Px(); - py = lv->Py(); - pz = lv->Pz(); - en = lv->Energy(); - - fastjet::PseudoJet input_part(px,py,pz,en); // create PseudoJet object - input_part.set_user_index(i); //label the particle into Fastjet algortihm - input_particles.push_back(input_part); // back of the input_particles vector - } // end loop - + if(fOpt==0) + { + TClonesArray *lvArray = fReader->GetMomentumArray(); + if(lvArray == 0) { cout << "Could not get the momentum array" << endl; return; } + Int_t nIn = lvArray->GetEntries(); + if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; } + fJets->SetNinput(nIn) ; // number of input objects + Float_t px,py,pz,en; + // load input vectors + for(Int_t i = 0; i < nIn; i++){ // loop for all input particles + TLorentzVector *lv = (TLorentzVector*) lvArray->At(i); + px = lv->Px(); + py = lv->Py(); + pz = lv->Pz(); + en = lv->Energy(); + + fastjet::PseudoJet input_part(px,py,pz,en); // create PseudoJet object + input_part.set_user_index(i); //label the particle into Fastjet algortihm + input_particles.push_back(input_part); // back of the input_particles vector + } // end loop + } + else { + TClonesArray* fUnit = fReader->GetUnitArray(); + if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; } + Int_t nCandidate = fReader->GetNumCandidate(); + Int_t nIn = fUnit->GetEntries(); + if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; } + fJets->SetNinput(nCandidate); // number of input objects // ME + // Information extracted from fUnitArray + // load input vectors and calculate total energy in array + Float_t pt,eta,phi,theta,px,py,pz,en; + Int_t ipart = 0; + for(Int_t i=0; iAt(i); + + if(uArray->GetUnitEnergy()>0.){ + + // It is not necessary anymore to cut on particle pt + pt = uArray->GetUnitEnergy(); + eta = uArray->GetUnitEta(); + phi = uArray->GetUnitPhi(); + theta = EtaToTheta(eta); + en = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta)); + px = TMath::Cos(phi)*pt; + py = TMath::Sin(phi)*pt; + pz = en*TMath::TanH(eta); + if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl; + + fastjet::PseudoJet input_part(px,py,pz,en); // create PseudoJet object + input_part.set_user_index(ipart); //label the particle into Fastjet algortihm + input_particles.push_back(input_part); // back of the input_particles vector + ipart++; + } + } // End loop on UnitArray + } // create an object that represents your choice of jet algorithm, and // the associated parameters @@ -116,8 +158,7 @@ void AliFastJetFinder::FindJets() fastjet::RecombinationScheme recomb_scheme = header->GetRecombScheme(); fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); fastjet::JetDefinition jet_def(algorithm, Rparam, recomb_scheme, strategy); - - + // create an object that specifies how we to define the area fastjet::AreaDefinition area_def; double ghost_etamax = header->GetGhostEtaMax(); @@ -130,71 +171,114 @@ void AliFastJetFinder::FindJets() fastjet::AreaType area_type = header->GetAreaType(); area_def = fastjet::AreaDefinition(area_type,ghost_spec); - // run the jet clustering with the above jet definition - fastjet::ClusterSequenceArea clust_seq(input_particles, jet_def, area_def); - - - // save a comment in the header - - TString comment = "Running FastJet algorithm with the following setup. "; - comment+= "Jet definition: "; - comment+= TString(jet_def.description()); - comment+= ". Area definition: "; - comment+= TString(area_def.description()); - comment+= ". Strategy adopted by FastJet: "; - comment+= TString(clust_seq.strategy_string()); - header->SetComment(comment); - if(debug){ - cout << "--------------------------------------------------------" << endl; - cout << comment << endl; - cout << "--------------------------------------------------------" << endl; - } - //header->PrintParameters(); - - - // extract the inclusive jets with pt > ptmin, sorted by pt - double ptmin = header->GetPtMin(); - vector inclusive_jets = clust_seq.inclusive_jets(ptmin); - - //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl; - - - //subtract background // =========================================== - // set the rapididty , phi range within which to study the background - double rap_max = header->GetRapMax(); - double rap_min = header->GetRapMin(); - double phi_max = header->GetPhiMax(); - double phi_min = header->GetPhiMin(); - fastjet::RangeDefinition range(rap_min, rap_max, phi_min, phi_max); - // subtract background - vector sub_jets = clust_seq.subtracted_jets(range,ptmin); - - // print out - //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n"; - //cout << "---------------------------------------\n"; - //cout << endl; - //printf(" ijet rap phi Pt area +- err\n"); - - // sort jets into increasing pt - vector jets = sorted_by_pt(sub_jets); - for (size_t j = 0; j < jets.size(); j++) { // loop for jets + if(bgMode) // BG subtraction + { + //***************************** JETS FINDING AND EXTRACTION + // run the jet clustering with the above jet definition + fastjet::ClusterSequenceArea clust_seq(input_particles, jet_def, area_def); -// double area = clust_seq.area(jets[j]); -// double area_error = clust_seq.area_error(jets[j]); - - //printf("%5u %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, area_error); + // save a comment in the header + + TString comment = "Running FastJet algorithm with the following setup. "; + comment+= "Jet definition: "; + comment+= TString(jet_def.description()); + comment+= ". Area definition: "; + comment+= TString(area_def.description()); + comment+= ". Strategy adopted by FastJet: "; + comment+= TString(clust_seq.strategy_string()); + header->SetComment(comment); + if(debug){ + cout << "--------------------------------------------------------" << endl; + cout << comment << endl; + cout << "--------------------------------------------------------" << endl; + } + //header->PrintParameters(); + + + // extract the inclusive jets with pt > ptmin, sorted by pt + double ptmin = header->GetPtMin(); + vector inclusive_jets = clust_seq.inclusive_jets(ptmin); + + //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl; + + + //subtract background // =========================================== + // set the rapididty , phi range within which to study the background + double rap_max = header->GetRapMax(); + double rap_min = header->GetRapMin(); + double phi_max = header->GetPhiMax(); + double phi_min = header->GetPhiMin(); + fastjet::RangeDefinition range(rap_min, rap_max, phi_min, phi_max); + + // subtract background + vector sub_jets = clust_seq.subtracted_jets(range,ptmin); + + // print out + //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n"; + //cout << "---------------------------------------\n"; + //cout << endl; + //printf(" ijet rap phi Pt area +- err\n"); + + // sort jets into increasing pt + vector jets = sorted_by_pt(sub_jets); + for (size_t j = 0; j < jets.size(); j++) { // loop for jets + + double area = clust_seq.area(jets[j]); + double area_error = clust_seq.area_error(jets[j]); + + printf("Jet found %5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n", (Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, area_error); // go to write AOD info - AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E()); - //cout << "Printing jet " << endl; - if(debug) aodjet.Print(""); - //cout << "Adding jet ... " ; - AddJet(aodjet); - //cout << "added \n" << endl; - - } // end loop for jets + AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E()); + //cout << "Printing jet " << endl; + if(debug) aodjet.Print(""); + //cout << "Adding jet ... " ; + AddJet(aodjet); + //cout << "added \n" << endl; + + } + } + else { // No BG subtraction + fastjet::ClusterSequence clust_seq(input_particles, jet_def); + + // save a comment in the header + + TString comment = "Running FastJet algorithm with the following setup. "; + comment+= "Jet definition: "; + comment+= TString(jet_def.description()); + comment+= ". Strategy adopted by FastJet: "; + comment+= TString(clust_seq.strategy_string()); + header->SetComment(comment); + if(debug){ + cout << "--------------------------------------------------------" << endl; + cout << comment << endl; + cout << "--------------------------------------------------------" << endl; + } + //header->PrintParameters(); + + // extract the inclusive jets with pt > ptmin, sorted by pt + double ptmin = header->GetPtMin(); + vector inclusive_jets = clust_seq.inclusive_jets(ptmin); + + //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl; + vector jets = sorted_by_pt(inclusive_jets); // Added by me + for (size_t j = 0; j < jets.size(); j++) { // loop for jets // Added by me + + printf("Jet found %5d %9.5f %8.5f %10.3f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp()); + + // go to write AOD info + AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E()); + //cout << "Printing jet " << endl; + if(debug) aodjet.Print(""); + //cout << "Adding jet ... " ; + AddJet(aodjet); + //cout << "added \n" << endl; + + } // end loop for jets + } + } //____________________________________________________________________________ @@ -282,7 +366,7 @@ void AliFastJetFinder::RunTest(const char* datafile) double area = clust_seq.area(jets[j]); double area_error = clust_seq.area_error(jets[j]); - printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",j,jets[j].rap(), + printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(), jets[j].phi(),jets[j].perp(), area, area_error); } cout << endl; @@ -300,3 +384,21 @@ void AliFastJetFinder::WriteJHeaderToFile() } //____________________________________________________________________________ + +Float_t AliFastJetFinder::EtaToTheta(Float_t arg) +{ + // return (180./TMath::Pi())*2.*atan(exp(-arg)); + return 2.*atan(exp(-arg)); + + +} + +//____________________________________________________________________________ + +void AliFastJetFinder::InitTask(TChain *tree) +{ + + printf("Fast jet finder initialization ******************"); + fReader->CreateTasks(tree); + +} diff --git a/JETAN/AliFastJetFinder.h b/JETAN/AliFastJetFinder.h index 55e60331192..a28ddb40c1c 100644 --- a/JETAN/AliFastJetFinder.h +++ b/JETAN/AliFastJetFinder.h @@ -30,7 +30,7 @@ #include #include "AliJetFinder.h" -#include "AliFastJetHeader.h" +#include "AliFastJetHeaderV1.h" using namespace std; @@ -42,12 +42,13 @@ class AliFastJetFinder : public AliJetFinder AliFastJetFinder(); ~AliFastJetFinder(); - void FindJets(); + void FindJets(); // others - void RunTest(const char* datafile); // a simple test - - void WriteJHeaderToFile(); - + void RunTest(const char* datafile); // a simple test + void WriteJHeaderToFile(); + Float_t EtaToTheta(Float_t arg); + void InitTask(TChain* tree); + protected: AliFastJetFinder(const AliFastJetFinder& rfj); AliFastJetFinder& operator = (const AliFastJetFinder& rsfj); diff --git a/JETAN/AliFastJetHeaderV1.cxx b/JETAN/AliFastJetHeaderV1.cxx new file mode 100644 index 00000000000..20a2ae12d92 --- /dev/null +++ b/JETAN/AliFastJetHeaderV1.cxx @@ -0,0 +1,89 @@ +/************************************************************************** + * 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. * + **************************************************************************/ + + +//--------------------------------------------------------------------- +// FastJet v2.3.4 finder algorithm interface +// Finder Header Class +// Author: Rafael.Diaz.Valdes@cern.ch +//--------------------------------------------------------------------- + +#include +#include + +#include "fastjet/ClusterSequenceArea.hh" +#include "fastjet/AreaDefinition.hh" +#include "fastjet/JetDefinition.hh" + +#include "AliFastJetHeaderV1.h" + +ClassImp(AliFastJetHeaderV1) + +//////////////////////////////////////////////////////////////////////// + +AliFastJetHeaderV1::AliFastJetHeaderV1(): + AliJetHeader("AliFastJetHeaderV1"), + fRparam(1.0), + fAlgorithm(fastjet::kt_algorithm), + fStrategy(fastjet::Best), + fRecombScheme(fastjet::BIpt_scheme), + fGhostEtaMax(2.0), + fGhostArea(0.05), + fActiveAreaRepeats(1), + fAreaType(fastjet::active_area), + fPtMin(5.0), + fPhiMin(0), + fPhiMax(TMath::TwoPi()), + fDebug(0), + fBGMode(0) +{ + // Constructor + + Double_t rapmax = fGhostEtaMax - fRparam; + Double_t rapmin = -fGhostEtaMax + fRparam; + SetRapRange(rapmin, rapmax); + +} + +//////////////////////////////////////////////////////////////////////// + +void AliFastJetHeaderV1::PrintParameters() const +{ + // prints out parameters of jet algorithm + + cout << "FastJet algorithm parameters:" << endl; + + cout << "-- Jet Definition --- " << endl; + cout << "R " << fRparam << endl; + cout << "Jet Algorithm " << fAlgorithm << endl; + cout << "Strategy " << fStrategy << endl; + cout << "Recombination Scheme " << fRecombScheme << endl; + + cout << "-- Ghosted Area Spec parameters --- " << endl; + cout << "Ghost Eta Max " << fGhostEtaMax << endl; + cout << "Ghost Area " << fGhostArea << endl; + cout << "Active Area Repeats " << fActiveAreaRepeats << endl; + + cout << "-- Area Definition parameters --- " << endl; + cout << "Area Type " << fAreaType << endl; + + cout << "-- Cluster Sequence Area parameters --- " << endl; + cout << "pt min " << fPtMin << endl; + + cout << "-- Range Definition parameters --- " << endl; + cout << " bkg rapidity range from " << fRapMin << " to " << fRapMax << endl; + cout << " bkg phi range from " << fPhiMin << " to " << fPhiMax << endl; + +} diff --git a/JETAN/AliFastJetHeaderV1.h b/JETAN/AliFastJetHeaderV1.h new file mode 100644 index 00000000000..a8e541d0e6a --- /dev/null +++ b/JETAN/AliFastJetHeaderV1.h @@ -0,0 +1,97 @@ +#ifndef ALIFASTJETHEADERV1_H +#define ALIFASTJETHEADERV1_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//--------------------------------------------------------------------- +// FastJet v2.3.4 finder algorithm interface +// Finder Header Class +// Author: Rafael.Diaz.Valdes@cern.ch +//--------------------------------------------------------------------- + + +#include "fastjet/ClusterSequenceArea.hh" +#include "fastjet/AreaDefinition.hh" +#include "fastjet/JetDefinition.hh" + +#include "AliJetHeader.h" + + +class AliFastJetHeaderV1 : public AliJetHeader +{ + public: + + AliFastJetHeaderV1(); + virtual ~AliFastJetHeaderV1() { } + + // Getters + Double_t GetRparam() const {return fRparam;} + fastjet::JetAlgorithm GetAlgorithm() const {return fAlgorithm;} + fastjet::Strategy GetStrategy() const {return fStrategy;} + fastjet::RecombinationScheme GetRecombScheme() const {return fRecombScheme;} + Double_t GetGhostEtaMax() const {return fGhostEtaMax;} + Double_t GetGhostArea() const {return fGhostArea;} + Int_t GetActiveAreaRepeats() const {return fActiveAreaRepeats;} + fastjet::AreaType GetAreaType() const {return fAreaType;} + Double_t GetPtMin() const {return fPtMin;} + Double_t GetRapMax() const {return fRapMax;} + Double_t GetRapMin() const {return fRapMin;} + Double_t GetPhiMax() const {return fPhiMax;} + Double_t GetPhiMin() const {return fPhiMin;} + // Added temporarily !!! To be removed if not necessary + Float_t GetMinCellEt() const {return fMinCellEt;} + Bool_t GetDebug() const {return fDebug;} + Bool_t GetBGMode() const {return fBGMode;} + + // Setters + void SetRparam(Double_t f) {fRparam = f;} + void SetAlgorithm(fastjet::JetAlgorithm f) {fAlgorithm = f;} + void SetStrategy(fastjet::Strategy f) {fStrategy = f;} + void SetRecombScheme(fastjet::RecombinationScheme f) {fRecombScheme = f;} + void SetGhostEtaMax(Double_t f) {fGhostEtaMax = f;} + void SetGhostArea(Double_t f) {fGhostArea = f;} + void SetActiveAreaRepeats(Int_t f) {fActiveAreaRepeats =f;} + void SetAreaType(fastjet::AreaType f) {fAreaType = f;} + void SetRapRange(Double_t fmin, Double_t fmax) {fRapMin = fmin; fRapMax = fmax;} + void SetPhiRange(Double_t fmin, Double_t fmax) {fPhiMin = fmin; fPhiMax = fmax;} + void SetPtMin(Double_t ptmin) {fPtMin = ptmin;} + void SetDebug(Bool_t debug) {fDebug = debug;} + void SetBGMode(Bool_t bgmode) {fBGMode = bgmode;} + + void SetComment(TString com) {fComment=com;} + void SetComment(const char* com) {AliJetHeader::SetComment(com);} + + // others + void PrintParameters() const; + + protected: + + //fastjet::JetDefinition parameters + Double_t fRparam; + fastjet::JetAlgorithm fAlgorithm; //fastjet::kt_algorithm + fastjet::Strategy fStrategy; //= fastjet::Best; + fastjet::RecombinationScheme fRecombScheme; // = fastjet::BIpt_scheme; + + //fastjet::GhostedAreaSpec parameters + Double_t fGhostEtaMax; + Double_t fGhostArea; + Int_t fActiveAreaRepeats; + + //fastjet::AreaDefinition parameters + fastjet::AreaType fAreaType; + + //fastjet::ClusterSequenceArea options parameters + Double_t fPtMin; //jets with pt > ptmin + Float_t fMinCellEt; // Min Et in one cell + + //fastjet::RangeDefinition parameters + Double_t fRapMax, fRapMin; // rapidity range of background sub + Double_t fPhiMax, fPhiMin; // phi range of background sub + Bool_t fDebug; // debug option + Bool_t fBGMode; // Do we subtract BG or not? + + ClassDef(AliFastJetHeaderV1,2) +}; + +#endif diff --git a/JETAN/AliJet.cxx b/JETAN/AliJet.cxx index d6940978b03..23e4e68739e 100644 --- a/JETAN/AliJet.cxx +++ b/JETAN/AliJet.cxx @@ -20,8 +20,10 @@ // Stores the output of a jet algorithm // Author: jgcn@mda.cinvestav.mx //--------------------------------------------------------------------- + #include +//#include #include #include @@ -39,7 +41,15 @@ AliJet::AliJet(): fJets(0), fEtaIn(0), fPhiIn(0), - fPtIn(0) + fPtIn(0), + fPtChPtCutIn(0), + fEnTotChPtCutIn(0), + fVectorSizeIn(0), + fDetIn(0), + fVPx(0), + fVPy(0), + fVPz(0) + // fVectorIn(0) { // Default constructor fJets = new TClonesArray("TLorentzVector",1000); @@ -50,6 +60,10 @@ AliJet::AliJet(): fPtFromSignal = TArrayF(); fMultiplicities = TArrayI(); fNCells = TArrayI(); + fPtChPtCutIn = TArrayF(); + fEnTotChPtCutIn = TArrayF(); + fVectorSizeIn = TArrayI(); + fDetIn = TArrayI(); } //////////////////////////////////////////////////////////////////////// @@ -244,6 +258,53 @@ void AliJet::SetPhiIn(Float_t* x) { if (fNInput>0) fPhiIn.Set(fNInput, x); } +//////////////////////////////////////////////////////////////////////// + +void AliJet::SetPtChargedPtCutIn(Float_t* x) +{ + if (fNInput>0) fPtChPtCutIn.Set(fNInput, x); +} +//////////////////////////////////////////////////////////////////////// + +void AliJet::SetEnTotChargedPtCutIn(Float_t* x) +{ + if (fNInput>0) fEnTotChPtCutIn.Set(fNInput, x); +} + +//////////////////////////////////////////////////////////////////////// + +void AliJet::SetDetectorFlagIn(Int_t* x) +{ + if (fNInput>0) fDetIn.Set(fNInput, x); +} + +//////////////////////////////////////////////////////////////////////// + +void AliJet::SetVectorSizeIn(Int_t* x) +{ + if (fNInput>0) fVectorSizeIn.Set(fNInput, x); +} + +//////////////////////////////////////////////////////////////////////// + +void AliJet::SetVectorPxIn(vector< vector > pxT) +{ + fVPx = pxT;; +} + +//////////////////////////////////////////////////////////////////////// + +void AliJet::SetVectorPyIn(vector< vector > pyT) +{ + fVPy = pyT;; +} + +//////////////////////////////////////////////////////////////////////// + +void AliJet::SetVectorPzIn(vector< vector > pzT) +{ + fVPz = pzT;; +} //////////////////////////////////////////////////////////////////////// @@ -285,6 +346,10 @@ void AliJet::ClearJets(Option_t *option) fEtaIn.Set(0); fPtIn.Set(0); fNCells.Set(0); + fPtChPtCutIn.Set(0); + fEnTotChPtCutIn.Set(0); + fVectorSizeIn.Set(0); + fDetIn.Set(0); } //////////////////////////////////////////////////////////////////////// diff --git a/JETAN/AliJet.h b/JETAN/AliJet.h index 84fb508b75c..74689e95edf 100644 --- a/JETAN/AliJet.h +++ b/JETAN/AliJet.h @@ -3,14 +3,16 @@ /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ - -/* $Id$ */ + //--------------------------------------------------------------------- // Jet class // Stores the output of a jet algorithm // Author: jgcn@mda.cinvestav.mx //--------------------------------------------------------------------- + +#include +#include #include #include @@ -27,21 +29,28 @@ class AliJet : public TObject ~AliJet(); // Getters - Int_t GetNinput() const {return fNInput;} - Int_t GetNJets() const {return fNJets;} + Int_t GetNinput() const { return fNInput; } + Int_t GetNJets() const {return fNJets;} TClonesArray* GetJets() const {return fJets;} - TArrayI GetInJet() const {return fInJet;} - TArrayI GetMultiplicities() const {return fMultiplicities;} - TArrayI GetNCells() const {return fNCells;} - TArrayF GetPtFromSignal() const {return fPtFromSignal;} - TArrayF GetEtaIn() const {return fEtaIn;} - TArrayF GetPhiIn() const {return fPhiIn;} - TArrayF GetPtIn() const {return fPtIn;} - Double_t GetEtAvg() const {return fEtAvg;} + TArrayI GetInJet() const {return fInJet;} + TArrayI GetMultiplicities() const {return fMultiplicities;} + TArrayI GetNCells() const {return fNCells;} + TArrayF GetPtFromSignal() const {return fPtFromSignal;} + TArrayF GetEtaIn() const { return fEtaIn; } + TArrayF GetPhiIn() const { return fPhiIn; } + TArrayF GetPtIn() const { return fPtIn; } + TArrayF GetPtChargedPtCutIn() const { return fPtChPtCutIn; } + TArrayF GetEnTotChargedPtCutIn() const {return fEnTotChPtCutIn; } + TArrayI GetVectorSizeIn() const { return fVectorSizeIn; } + vector< vector > GetVectorPxIn() const { return fVPx; } + vector< vector > GetVectorPyIn() const { return fVPy; } + vector< vector > GetVectorPzIn() const { return fVPz; } + TArrayI GetDetectorFlagIn() const { return fDetIn; } + Double_t GetEtAvg() const { return fEtAvg; } TLorentzVector* GetJet(Int_t i); - Int_t GetMultiplicity(Int_t i) const; - Int_t GetNCell(Int_t i) const; + Int_t GetMultiplicity(Int_t i) const; + Int_t GetNCell(Int_t i) const; Double_t GetPx(Int_t i); Double_t GetPy(Int_t i); Double_t GetPz(Int_t i); @@ -63,8 +72,14 @@ class AliJet : public TObject void SetPhiIn(Float_t* phi); void SetPtIn(Float_t* pt); void SetInJet(Int_t* idx); + void SetPtChargedPtCutIn(Float_t* pt2T); + void SetEnTotChargedPtCutIn(Float_t* en2T); + void SetVectorSizeIn(Int_t* vectT); + void SetVectorPxIn(vector< vector > pxT); + void SetVectorPyIn(vector< vector > pyT); + void SetVectorPzIn(vector< vector > pzT); + void SetDetectorFlagIn(Int_t* detT); void SetEtAvg(Double_t et) { fEtAvg = et; } - // others Bool_t OutOfRange(Int_t i, const char *s) const; void ClearJets(Option_t *option=""); @@ -75,21 +90,28 @@ class AliJet : public TObject AliJet(const AliJet& rJet); AliJet& operator = (const AliJet& rhs); - Int_t fNInput; // number of input objects - Int_t fNJets; // number of jets found - Double_t fEtAvg; // average background et per cell + Int_t fNInput; // number of input objects + Int_t fNJets; // number of jets found + Double_t fEtAvg; // average background et per cell + + TArrayI fInJet; // i-input object belongs to k-jet + TArrayI fMultiplicities; // Multiplicity of each jet + TArrayI fNCells; // Number of cells in jet + TArrayF fPtFromSignal; // percentage of pt from signal + TClonesArray* fJets; // 4-momenta of jets - TArrayI fInJet; // i-input object belongs to k-jet - TArrayI fMultiplicities; // Multiplicity of each jet - TArrayI fNCells; // Number of cells in jet - TArrayF fPtFromSignal; // percentage of pt from signal - TClonesArray* fJets; // 4-momenta of jets + TArrayF fEtaIn; // Arrays of input particles kine:Eta + TArrayF fPhiIn; // Arrays of input particles kine:Phi + TArrayF fPtIn; // Arrays of input particles kine:Pt + TArrayF fPtChPtCutIn; // Arrays of input particles kin:Pt Charged with pt cut + TArrayF fEnTotChPtCutIn; // Arrays of total energy with pt cut on charged + cut min on cell + TArrayI fVectorSizeIn; // Arrays of number of charged tracks in each unitArray + TArrayI fDetIn; // Arrays of detector type of each UnitArray + vector< vector > fVPx; //|| + vector< vector > fVPy; //|| + vector< vector > fVPz; //|| - TArrayF fEtaIn; // Arrays of input particles kine:Eta - TArrayF fPhiIn; // Arrays of input particles kine:Phi - TArrayF fPtIn; // Arrays of input particles kine:Pt - - ClassDef(AliJet,1) + ClassDef(AliJet,2) }; #endif diff --git a/JETAN/AliJetESDReader.cxx b/JETAN/AliJetESDReader.cxx index 1ac583d8f6a..055f1e7467c 100755 --- a/JETAN/AliJetESDReader.cxx +++ b/JETAN/AliJetESDReader.cxx @@ -17,7 +17,7 @@ // Jet ESD Reader // ESD reader for jet analysis // Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch) -// Magali Estienne +// Magali Estienne //------------------------------------------------------------------------- // --- Standard library --- @@ -28,13 +28,15 @@ #include #include #include -#include +#include "TTask.h" +#include "TTree.h" +#include "TFile.h" #include #include #include #include -#include - +#include +#include // --- AliRoot header files --- #include "AliJetESDReader.h" @@ -46,6 +48,7 @@ #include "AliJetFillUnitArrayTracks.h" #include "AliJetFillUnitArrayEMCalDigits.h" #include "AliJetUnitArray.h" +#include "AliAnalysisTask.h" ClassImp(AliJetESDReader) @@ -53,6 +56,7 @@ AliJetESDReader::AliJetESDReader(): AliJetReader(), fGeom(0), fChain(0x0), + fTree(0x0), fESD(0x0), fHadCorr(0x0), fTpcGrid(0x0), @@ -64,6 +68,8 @@ AliJetESDReader::AliJetESDReader(): fGrid4(0), fPtCut(0), fHCorrection(0), + fECorrection(0), + fEFlag(kFALSE), fNumUnits(0), fDebug(0), fMass(0), @@ -73,7 +79,9 @@ AliJetESDReader::AliJetESDReader(): fDZ(0), fNeta(0), fNphi(0), - fArrayInitialised(0) + fArrayInitialised(0), + fRefArray(0x0), + fProcId(kFALSE) { // Constructor } @@ -83,6 +91,7 @@ AliJetESDReader::~AliJetESDReader() { // Destructor delete fChain; + delete fTree; delete fESD; delete fTpcGrid; delete fEmcalGrid; @@ -117,7 +126,6 @@ void AliJetESDReader::OpenInputFiles() if (strstr(name,pattern)){ printf("Adding %s\n",name); char path[256]; - // sprintf(path,"%s/%s/AliESDs.root",dirName,name); sprintf(path,"%s/%s/AliESDs.root",dirName,name); fChain->AddFile(path); a++; @@ -132,7 +140,7 @@ void AliJetESDReader::OpenInputFiles() int nMax = fChain->GetEntries(); printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax); - + // set number of events in header if (fReaderHeader->GetLastEvent() == -1) fReaderHeader->SetLastEvent(nMax); @@ -140,6 +148,7 @@ void AliJetESDReader::OpenInputFiles() Int_t nUsr = fReaderHeader->GetLastEvent(); fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr)); } + } //____________________________________________________________________________ @@ -170,7 +179,7 @@ Bool_t AliJetESDReader::FillMomentumArray() // get number of tracks in event (for the loop) nt = fESD->GetNumberOfTracks(); printf("Fill Momentum Array %5d \n", nt); - + // temporary storage of signal and pt cut flag Int_t* sflag = new Int_t[nt]; Int_t* cflag = new Int_t[nt]; @@ -180,7 +189,7 @@ Bool_t AliJetESDReader::FillMomentumArray() Float_t etaMin = fReaderHeader->GetFiducialEtaMin(); Float_t etaMax = fReaderHeader->GetFiducialEtaMax(); - //loop over tracks + //loop over tracks in ESD for (Int_t it = 0; it < nt; it++) { AliESDtrack *track = fESD->GetTrack(it); UInt_t status = track->GetStatus(); @@ -205,6 +214,7 @@ Bool_t AliJetESDReader::FillMomentumArray() if (pt > ptMin) cflag[goodTrack]=1; // pt cut goodTrack++; } + // set the signal flags fSignalFlag.Set(goodTrack,sflag); fCutFlag.Set(goodTrack,cflag); @@ -213,28 +223,37 @@ Bool_t AliJetESDReader::FillMomentumArray() delete[] cflag; return kTRUE; + } //____________________________________________________________________________ -void AliJetESDReader::CreateTasks() +void AliJetESDReader::CreateTasks(TChain* tree) { + // + // For reader task initialization + // + fDebug = fReaderHeader->GetDebug(); fDZ = fReaderHeader->GetDZ(); + fTree = tree; - // Init EMCAL geometry and create UnitArray object + // Init EMCAL geometry SetEMCALGeometry(); + // Init parameters InitParameters(); + // Create and init unit array InitUnitArray(); + // Create global reader task for analysis fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder"); + // Create a task for to fill the charged particle information fFillUAFromTracks = new AliJetFillUnitArrayTracks(); fFillUAFromTracks->SetReaderHeader(fReaderHeader); fFillUAFromTracks->SetGeom(fGeom); fFillUAFromTracks->SetTPCGrid(fTpcGrid); fFillUAFromTracks->SetEMCalGrid(fEmcalGrid); - if(fDZ) - { + { // Calo dead zones inclusion fFillUAFromTracks->SetGrid0(fGrid0); fFillUAFromTracks->SetGrid1(fGrid1); fFillUAFromTracks->SetGrid2(fGrid2); @@ -243,12 +262,15 @@ void AliJetESDReader::CreateTasks() } fFillUAFromTracks->SetHadCorrection(fHCorrection); fFillUAFromTracks->SetHadCorrector(fHadCorr); + // Create a task for to fill the neutral particle information fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits(); fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader); fFillUAFromEMCalDigits->SetGeom(fGeom); fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid); fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid); - // fFillUnitArray->Add(fFillUAFromTracks); + fFillUAFromEMCalDigits->SetEleCorrection(fECorrection); + // Add the task to global task + fFillUnitArray->Add(fFillUAFromTracks); fFillUnitArray->Add(fFillUAFromEMCalDigits); fFillUAFromTracks->SetActive(kFALSE); fFillUAFromEMCalDigits->SetActive(kFALSE); @@ -259,49 +281,35 @@ void AliJetESDReader::CreateTasks() } //____________________________________________________________________________ -//void AliJetESDReader::ExecTasks(Int_t event) -Bool_t AliJetESDReader::ExecTasks(Int_t /*event*/) +Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray) { - // clear momentum array - Int_t nEntRef = fRefArray->GetEntries(); - - for(Int_t i=0; iAt(i))->SetUnitTrackID(0); - ((AliJetUnitArray*)fRefArray->At(i))->SetUnitEnergy(0.); - ((AliJetUnitArray*)fRefArray->At(i))->SetUnitCutFlag(kPtSmaller); - ((AliJetUnitArray*)fRefArray->At(i))->SetUnitDetectorFlag(kTpc); - ((AliJetUnitArray*)fRefArray->At(i))->SetUnitFlag(kOutJet); - } + // + // Reader task execussion + // + + fProcId = procid; + fRefArray = refArray; + vector vtmp(3); + // clear momentum array ClearArray(); fDebug = fReaderHeader->GetDebug(); fOpt = fReaderHeader->GetDetector(); - // InitParameters(); if(!fESD) { return kFALSE; } - /* - // get event from chain - // For TSelectors - // fChain->GetTree()->GetEntry(event); - // For interactive process - // fChain->GetEntry(event); - fChain->GetEvent(event); - */ - // TPC only or Digits+TPC or Clusters+TPC if(fOpt%2==!0 && fOpt!=0){ fFillUAFromTracks->SetESD(fESD); fFillUAFromTracks->SetActive(kTRUE); fFillUAFromTracks->SetUnitArray(fUnitArray); fFillUAFromTracks->SetRefArray(fRefArray); - // fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed !!! - // Incompatibility with Andrei's analysis framework - fFillUAFromTracks->Exec("tpc"); + fFillUAFromTracks->SetProcId(fProcId); + // fFillUAFromTracks->ExecuteTask("tpc"); // Temporarily changed + fFillUAFromTracks->Exec("tpc"); // Temporarily added if(fOpt==1){ fNumCandidate = fFillUAFromTracks->GetMult(); fNumCandidateCut = fFillUAFromTracks->GetMultCut(); @@ -314,48 +322,71 @@ Bool_t AliJetESDReader::ExecTasks(Int_t /*event*/) fFillUAFromEMCalDigits->SetActive(kTRUE); fFillUAFromEMCalDigits->SetUnitArray(fUnitArray); fFillUAFromEMCalDigits->SetRefArray(fRefArray); + fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId()); fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult()); fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut()); - fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily changed !!! + fFillUAFromEMCalDigits->Exec("digits"); // Temporarily added fNumCandidate = fFillUAFromEMCalDigits->GetMult(); fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut(); } - // fFillUnitArray->ExecuteTask(); // => Temporarily commented + // fFillUnitArray->ExecuteTask(); // Temporarily commented return kTRUE; } //____________________________________________________________________________ -void AliJetESDReader::SetEMCALGeometry() +Bool_t AliJetESDReader::SetEMCALGeometry() { - // Define EMCAL geometry to be able to read ESDs - fGeom = AliJetDummyGeo::GetInstance(); - if (fGeom == 0) - fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL"); + // + // Set the EMCal Geometry + // - // To be setted to run some AliEMCALGeometry functions - TGeoManager::Import("geometry.root"); - // fGeom->GetTransformationForSM(); - printf("\n EMCal Geometry set ! \n"); + if (!fTree->GetFile()) + return kFALSE; -} + TString geomFile(fTree->GetFile()->GetName()); + geomFile.ReplaceAll("AliESDs", "geometry"); + + // temporary workaround for PROOF bug #18505 + geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root"); + if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data()); + // Define EMCAL geometry to be able to read ESDs + fGeom = AliJetDummyGeo::GetInstance(); + if (fGeom == 0) + fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL"); + + // To be setted to run some AliEMCALGeometry functions + TGeoManager::Import(geomFile); + fGeom->GetTransformationForSM(); + printf("\n EMCal Geometry set ! \n"); + + return kTRUE; +} //____________________________________________________________________________ void AliJetESDReader::InitParameters() { - // Initialise parameters - fHCorrection = 0; // For hadron correction - fHadCorr = 0; // For hadron correction - fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL - if(fDebug>1) printf("\n EMCal parameters initiated ! \n"); + // Initialise parameters + fOpt = fReaderHeader->GetDetector(); + fHadCorr = 0; // For hadron correction + if(fEFlag==kFALSE){ + if(fOpt==0 || fOpt==1) + fECorrection = 0; // For electron correction + else fECorrection = 1; // For electron correction + } + fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL + if(fDebug>1) printf("\n EMCal parameters initiated ! \n"); } //____________________________________________________________________________ void AliJetESDReader::InitUnitArray() { - //Initialises unit arrays + // + // Create and Initialise unit arrays + // + Int_t nElements = fTpcGrid->GetNEntries(); Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.; if(fArrayInitialised) fUnitArray->Delete(); @@ -366,12 +397,11 @@ void AliJetESDReader::InitUnitArray() // detector flag, in/out jet, pt cut, mass, cluster ID) for(Int_t nBin = 1; nBin < nElements+1; nBin++) { - // fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi); fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta); phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi); Deta = fTpcGrid->GetDeta(); Dphi = fTpcGrid->GetDphi(); - new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); } } @@ -438,7 +468,7 @@ void AliJetESDReader::InitUnitArray() phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi); Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values - new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); } else { if(nBin>=fNumUnits && nBinGetDeta(); Dphi = fTpcGrid->GetDphi(); - new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); } else { if(fDZ) { if(nBin>=fNumUnits+nElements && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta); Deta = fGrid0->GetDeta(); Dphi = fGrid0->GetDphi(); - new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); } else if(nBin>=fNumUnits+nElements+n0 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta); Deta = fGrid1->GetDeta(); Dphi = fGrid1->GetDphi(); - new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); } else if(nBin>=fNumUnits+nElements+n0+n1 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta); Deta = fGrid2->GetDeta(); Dphi = fGrid2->GetDphi(); - new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); } else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta); Deta = fGrid3->GetDeta(); Dphi = fGrid3->GetDphi(); - new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); } else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBinGetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta); Deta = fGrid4->GetDeta(); Dphi = fGrid4->GetDphi(); - new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1); + new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1); } } } // end if(fDZ) diff --git a/JETAN/AliJetESDReader.h b/JETAN/AliJetESDReader.h index 1d14954d61c..d42595c8718 100755 --- a/JETAN/AliJetESDReader.h +++ b/JETAN/AliJetESDReader.h @@ -11,13 +11,18 @@ //========================================================================= // Modified in order to use a fUnitArray object instead of a fMomentumArray // Includes EMCal Geometry, fUnitArray, grid objects and tools for Hadron correction -// Author : magali.estienne@ires.in2p3.fr +// Author : magali.estienne@subatech.in2p3.fr //--------------------------------------------------------------------- +#include + #include "AliJetReader.h" #include "AliJetUnitArray.h" #include "AliJetGrid.h" + +class TRefArray; class AliJetESDReaderHeader; +class AliEMCALGeometry; class AliJetDummyGeo; class AliJetHadronCorrection; class AliJetUnitArray; @@ -30,32 +35,29 @@ class AliJetESDReader : public AliJetReader AliJetESDReader(); virtual ~AliJetESDReader(); - // Getters - Float_t GetTrackMass() const {return fMass;} // returns mass of the track - Int_t GetTrackSign() const {return fSign;} // returns sign of the track + Bool_t FillMomentumArray(); + void OpenInputFiles(); + void InitUnitArray(); + void CreateTasks(TChain* tree); + Bool_t ExecTasks(Bool_t procid, TRefArray* refArray); + // Getters + Float_t GetTrackMass() const {return fMass;} // returns mass of the track + Int_t GetTrackSign() const {return fSign;} // returns sign of the track + // Setters - Bool_t FillMomentumArray(); - void OpenInputFiles(); - void InitUnitArray(); - void CreateTasks(); - // void ExecTasks(Int_t event); - Bool_t ExecTasks(Int_t event); - void SetInputEvent(TObject* esd, TObject* aod, TObject* mc); - virtual void SetTPCGrid(AliJetGrid *grid) {fTpcGrid = grid;} - virtual void SetEMCalGrid(AliJetGrid *grid) {fEmcalGrid = grid;} + void SetInputEvent(TObject* esd, TObject* aod, TObject* mc); + void SetTPCGrid(AliJetGrid *grid) {fTpcGrid = grid;} + void SetEMCalGrid(AliJetGrid *grid) {fEmcalGrid = grid;} // Correction of hadronic energy - virtual void SetHadronCorrection(Int_t flag = 1) {fHCorrection = flag;} - virtual void SetHadronCorrector(AliJetHadronCorrectionv1* corr) {fHadCorr = corr;} - private: - void SetEMCALGeometry(); - void InitParameters(); - AliJetESDReader(const AliJetESDReader &det); - AliJetESDReader &operator=(const AliJetESDReader &det); + void SetHadronCorrection(Int_t flag = 1) {fHCorrection = flag;} + void SetHadronCorrector(AliJetHadronCorrectionv1* corr) {fHadCorr = corr;} + void SetElectronCorrection(Int_t flag = 1) {fECorrection = flag; fEFlag=kTRUE;} protected: AliJetDummyGeo *fGeom; //! EMCAL Geometry TChain *fChain; //! chain for reconstructed tracks + TChain *fTree; //! tree for reconstructed tracks AliESDEvent *fESD; //! pointer to esd AliJetHadronCorrectionv1 *fHadCorr; //! Pointer to Hadron Correction Object AliJetGrid *fTpcGrid; //! Pointer to grid object @@ -67,18 +69,30 @@ class AliJetESDReader : public AliJetReader AliJetGrid *fGrid4; // Pointer to grid object Float_t fPtCut; // Pt cut for tracks to minimise background contribution Int_t fHCorrection; // Hadron correction flag + Int_t fECorrection; // Electron correction flag + Bool_t fEFlag; // if (fEFlag == kFALSE) => fECorrection automatically setted Int_t fNumUnits; // Number of units in the unit object array // (same as num towers in EMCAL) - Int_t fDebug; // Debug option - Float_t fMass; // Particle mass - Int_t fSign; // Particle sign + Int_t fDebug; //! Debug option + Float_t fMass; // Particle mass + Int_t fSign; // Particle sign Int_t fNIn; // Number of Array filled in UnitArray Int_t fOpt; // Detector to be used for jet reconstruction Bool_t fDZ; // Use or not dead zones Int_t fNeta; // Number of bins in eta of tpc grid Int_t fNphi; // Number of bins in phi of tpc grid Bool_t fArrayInitialised; // To check that array of units is initialised - + TRefArray *fRefArray; // array of digit position and energy + Bool_t fProcId; // Bool_t for TProcessID synchronization + + private: + Bool_t SetEMCALGeometry(); + void InitParameters(); + AliJetESDReader(const AliJetESDReader &det); + AliJetESDReader &operator=(const AliJetESDReader &det); + + + ClassDef(AliJetESDReader,1) }; diff --git a/JETAN/AliJetFillUnitArrayEMCalDigits.cxx b/JETAN/AliJetFillUnitArrayEMCalDigits.cxx index 303d9b001f0..88dba8965ef 100644 --- a/JETAN/AliJetFillUnitArrayEMCalDigits.cxx +++ b/JETAN/AliJetFillUnitArrayEMCalDigits.cxx @@ -13,15 +13,14 @@ * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ - -//------------------------------------------------------------- -// Fill Unit Array with EMCal information + +// Fill Unit Array // Called by ESD reader for jet analysis // Author: Magali Estienne (magali.estienne@ires.in2p3.fr) -//------------------------------------------------------------- // --- Standard library --- #include +#include // --- ROOT system --- #include @@ -39,19 +38,27 @@ #include "AliJetReader.h" #include "AliJetESDReader.h" #include "AliJetESDReaderHeader.h" -//#include "AliESD.h" #include "AliESDEvent.h" +#include "AliESDVertex.h" #include "AliJetDummyGeo.h" #include "AliESDCaloCluster.h" +#include "AliESDCaloCells.h" #include "AliJetUnitArray.h" #include "AliJetFillUnitArrayEMCalDigits.h" +// Remove CDB dependence under construction +//#include "AliEMCALCalibData.h" +//#include "AliCDBManager.h" + +//class AliCDBStorage; +//#include "AliCDBEntry.h" + ClassImp(AliJetFillUnitArrayEMCalDigits) //_____________________________________________________________________________ AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits() : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"), - fESD(0), + fESD(0), fNumUnits(0), fEtaMinCal(0.), fEtaMaxCal(0.), @@ -59,28 +66,34 @@ AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits() fPhiMaxCal(0.), fNIn(0), fOpt(0), + fCluster(0), fDebug(0), fNCEMCAL(0), fNCPHOS(0), fNCCalo(0), fTPCGrid(0x0), fEMCalGrid(0x0), + fECorrection(0), fReaderHeader(0x0), fMomentumArray(0x0), fUnitArray(0x0), fRefArray(0x0), + fProcId(kFALSE), fGeom(0x0), fClus(0x0), fNDigitEmcal(0), - fNDigitEmcalCut(0) + fNDigitEmcalCut(0), + fCalibData(0x0), + fADCchannelECA(0), + fADCpedestalECA(0) { // constructor } //_____________________________________________________________________________ -AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent */*esd*/) +AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent *esd) : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"), - fESD(0), + fESD(esd), fNumUnits(0), fEtaMinCal(0.), fEtaMaxCal(0.), @@ -88,24 +101,102 @@ AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent */*es fPhiMaxCal(0.), fNIn(0), fOpt(0), + fCluster(0), fDebug(0), fNCEMCAL(0), fNCPHOS(0), fNCCalo(0), fTPCGrid(0x0), fEMCalGrid(0x0), + fECorrection(0), fReaderHeader(0x0), fMomentumArray(0x0), fUnitArray(0x0), fRefArray(0x0), + fProcId(kFALSE), fGeom(0x0), fClus(0x0), fNDigitEmcal(0), - fNDigitEmcalCut(0) + fNDigitEmcalCut(0), + fCalibData(0x0), + fADCchannelECA(0), + fADCpedestalECA(0) { - // constructor + // constructor 2 +} + +//_____________________________________________________________________________ +AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(const AliJetFillUnitArrayEMCalDigits &det) + : TTask(det),//"AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"), + fESD(det.fESD), + fNumUnits(det.fNumUnits), + fEtaMinCal(det.fEtaMinCal), + fEtaMaxCal(det.fEtaMaxCal), + fPhiMinCal(det.fPhiMinCal), + fPhiMaxCal(det.fPhiMaxCal), + fNIn(det.fNIn), + fOpt(det.fOpt), + fCluster(det.fCluster), + fDebug(det.fDebug), + fNCEMCAL(det.fNCEMCAL), + fNCPHOS(det.fNCPHOS), + fNCCalo(det.fNCCalo), + fTPCGrid(det.fTPCGrid), + fEMCalGrid(det.fEMCalGrid), + fECorrection(det.fECorrection), + fReaderHeader(det.fReaderHeader), + fMomentumArray(det.fMomentumArray), + fUnitArray(det.fUnitArray), + fRefArray(det.fRefArray), + fProcId(det.fProcId), + fGeom(det.fGeom), + fClus(det.fClus), + fNDigitEmcal(det.fNDigitEmcal), + fNDigitEmcalCut(det.fNDigitEmcalCut), + fCalibData(det.fCalibData), + fADCchannelECA(det.fADCchannelECA), + fADCpedestalECA(det.fADCpedestalECA) +{ + // Copy constructor } +//_____________________________________________________________________________ +AliJetFillUnitArrayEMCalDigits& AliJetFillUnitArrayEMCalDigits::operator=(const AliJetFillUnitArrayEMCalDigits& other) +{ + // Assignment + + fESD = other.fESD; + fNumUnits = other.fNumUnits; + fEtaMinCal = other.fEtaMinCal; + fEtaMaxCal = other.fEtaMaxCal; + fPhiMinCal = other.fPhiMinCal; + fPhiMaxCal = other.fPhiMaxCal; + fNIn = other.fNIn; + fOpt = other.fOpt; + fCluster = other.fCluster; + fDebug = other.fDebug; + fNCEMCAL = other.fNCEMCAL; + fNCPHOS = other.fNCPHOS; + fNCCalo = other.fNCCalo; + fTPCGrid = other.fTPCGrid; + fEMCalGrid = other.fEMCalGrid; + fECorrection = other.fECorrection; + fReaderHeader = other.fReaderHeader; + fMomentumArray = other.fMomentumArray; + fUnitArray = other.fUnitArray; + fRefArray = other.fRefArray; + fProcId = other.fProcId; + fGeom = other.fGeom; + fClus = other.fClus; + fNDigitEmcal = other.fNDigitEmcal; + fNDigitEmcalCut = other.fNDigitEmcalCut; + fCalibData = other.fCalibData; + fADCchannelECA = other.fADCchannelECA; + fADCpedestalECA = other.fADCpedestalECA; + + return (*this); + +} //____________________________________________________________________________ void AliJetFillUnitArrayEMCalDigits::InitParameters() @@ -118,6 +209,10 @@ void AliJetFillUnitArrayEMCalDigits::InitParameters() fPhiMaxCal = fGeom->GetArm1PhiMax(); fClus = 0; + // Get calibration parameters from file or digitizer default values. + // Under construction + // GetCalibrationParameters() ; + if(fDebug>1) printf("\n EMCAL parameters initiated ! \n"); } @@ -133,108 +228,259 @@ void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /*option*/) { // // Main method. - // Explain - - fDebug = fReaderHeader->GetDebug(); - fOpt = fReaderHeader->GetDetector(); - + // + + fDebug = fReaderHeader->GetDebug(); + fOpt = fReaderHeader->GetDetector(); + fCluster = fReaderHeader->GetCluster(); + // Init parameters InitParameters(); - - // Get number of clusters from EMCAL - Int_t nDigitTot = 0; Int_t goodDigit = 0; - Int_t beg = 0; - Int_t end = 0; - Float_t ptMin = fReaderHeader->GetPtCut(); - - // Loop over calo clusters - //------------------------------------------------------------------ - Int_t type = 0; - Int_t index = 0; + Int_t index = 0; + + if(!fCluster) { // Keep all digit information + // Loop over all cell information + //------------------------------------------------------------------ + AliESDCaloCells &cells = *(fESD->GetEMCALCells()); + Int_t ncell = cells.GetNumberOfCells() ; + + for (Int_t icell= 0; icell < ncell; icell++) { + Int_t digitID = cells.GetCellNumber(icell); + Double_t digitAmp = cells.GetAmplitude(icell); + Float_t digitEn = digitAmp*0.0153; // Last correct + // Float_t digitEn = Calibrate((Int_t)digitAmp,digitID); + + Float_t etaD=-10., phiD=-10.; + fGeom->EtaPhiFromIndex(digitID,etaD,phiD); + // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD); - // Total number of EMCAL cluster - end = fESD->GetNumberOfCaloClusters(); + phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD); + + Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); - for(Int_t j = beg; j < end; j++) { - fClus = fESD->GetCaloCluster(j); - if(!fClus->IsEMCAL()) continue; + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID); - type = fClus->GetClusterType(); - index = fClus->GetID(); - nDigitTot = fClus->GetNumberOfDigits(); + if(uArray->GetUnitEnergy() == 0.) goodDigit++; + uArray->SetUnitTrackID(digitID); + + Float_t unitEnergy = 0.; + Bool_t ok = kFALSE; + unitEnergy = uArray->GetUnitEnergy(); + + if(unitEnergy==0){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + fNDigitEmcal++; + ok = kTRUE; + } + + // Detector flag + if(unitEnergy>0.) + uArray->SetUnitDetectorFlag(kAll); + else uArray->SetUnitDetectorFlag(kEmcal); - // Keep clusters or pseudo clusters - if (type != AliESDCaloCluster::kEMCALClusterv1) continue; - // if (type != AliESDCaloCluster::kPseudoCluster) continue; + uArray->SetUnitEnergy(unitEnergy+digitEt); - // Get the digit index and the digit information - //============================================================ + uArray->SetUnitCutFlag(kPtHigher); + + // To be modified !!! + uArray->SetUnitSignalFlag(kGood); + + // This is for jet multiplicity + uArray->SetUnitClusterID(index); + + if(fDebug > 1) printf("goodDigit : %d\n", goodDigit); + + } // End loop over cells + + } // end if !fCluster + else { // Keep digit information from clusterization + + // Loop over calo clusters + //------------------------------------------------------------------ + + //select EMCAL clusters only + TRefArray * caloClusters = new TRefArray(); + fESD->GetEMCALClusters(caloClusters); + + // Total number of EMCAL cluster + Int_t nclus = caloClusters->GetEntries() ; + Int_t beg = 0; + Float_t pos[3] ; + + // Get reconstructed vertex position + Double_t vertex_position[3] ; + fESD->GetVertex()->GetXYZ(vertex_position) ; + + // Get CaloCells + AliESDCaloCells &cells= *(fESD->GetEMCALCells()); + + for(Int_t j = beg; j < nclus; j++) { // loop over clusters + // Retrieve cluster from esd + AliESDCaloCluster *fClus = (AliESDCaloCluster *) caloClusters->At(j) ; + + // Get the cluster info + Float_t energy = fClus->E() ; + Int_t iprim = fClus->GetLabel(); + Int_t trackIndex = fClus->GetTrackMatched(); + + fClus->GetPosition(pos) ; + TVector3 vpos(pos[0],pos[1],pos[2]) ; + TLorentzVector p ; + fClus->GetMomentum(p,vertex_position); + + Int_t digMult = fClus->GetNCells() ; + UShort_t *digID = fClus->GetCellsAbsId() ; + //Print cluster info + if(fDebug>2) cout<<"Cluster "<< j <<"; digits mult "<GetClusterType() + <<"; Energy "<GetCellNumber(i) ; + Double_t digitAmp = cells.GetCellAmplitude(digitID) ; - // Get number of digits in a cluster - Int_t nD = fClus->GetNumberOfDigits(); + // Calibration for an energy in GeV + Float_t digitEn = digitAmp*0.0153; + // Float_t digitEn = Calibrate((Int_t)digitAmp,digitID); - TArrayS *digID = fClus->GetDigitIndex(); - TArrayS *digEnergy = fClus->GetDigitAmplitude(); - Float_t *digitEnergy = new Float_t[nD]; - // Float_t digitEn = 0.; + Float_t etaD=-10., phiD=-10.; + fGeom->EtaPhiFromIndex(digitID,etaD,phiD); + // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD); - // Loop over digits - for(Int_t k=0; kAt(k); - // Calibration for an energy in GeV - digitEnergy[k] = (Float_t)digEnergy->At(k)/500.; - - // Second method to extract eta, phi positions of a digit - //================================================================= - - Float_t etaD=-10., phiD=-10.; - fGeom->EtaPhiFromIndex(idF,etaD,phiD); - // fEMCalGrid->GetEtaPhiFromIndex2(idF,phiD,etaD); - phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD); - - Float_t etDigit = digitEnergy[k]*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); - - AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idF); - if(uArray->GetUnitEnergy() == 0.) goodDigit++; - - Float_t unitEnergy = 0.; - Bool_t ok = kFALSE; - unitEnergy = uArray->GetUnitEnergy(); - if(unitEnergy==0){ - fRefArray->Add(uArray); - fNDigitEmcal++; - ok = kTRUE; - } - uArray->SetUnitEnergy(unitEnergy+etDigit); - // Put a cut flag - if(uArray->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); - else { + phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD); + + Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD))); + + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID); + if(uArray->GetUnitEnergy() == 0.) goodDigit++; + uArray->SetUnitTrackID(digitID); + + Float_t unitEnergy = 0.; + Bool_t ok = kFALSE; + unitEnergy = uArray->GetUnitEnergy(); + + if(unitEnergy==0){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + fNDigitEmcal++; + ok = kTRUE; + } + + // Detector flag + if(unitEnergy>0.) + uArray->SetUnitDetectorFlag(kAll); + else uArray->SetUnitDetectorFlag(kEmcal); + + uArray->SetUnitEnergy(unitEnergy+digitEt); uArray->SetUnitCutFlag(kPtHigher); - if(ok) fNDigitEmcalCut++; - } - // Detector flag - if(unitEnergy>0.) - uArray->SetUnitDetectorFlag(kAll); - else uArray->SetUnitDetectorFlag(kEmcal); - - // This is for jet multiplicity - uArray->SetUnitClusterID(index); - - if(fDebug > 12) printf("goodDigit : %d\n", goodDigit); - - } // End loop over digits - - } // End loop over clusters - + + // To be modified !!! + uArray->SetUnitSignalFlag(kGood); + + // This is for jet multiplicity + uArray->SetUnitClusterID(index); + + } // End loop over cells + } // End loop over clusters + } // end else + fNIn += goodDigit; + if(fDebug>1) + { + printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"); + printf("goodDigit : %d\n", goodDigit); + } + } +// //____________________________________________________________________________ +// void AliJetFillUnitArrayEMCalDigits::GetCalibrationParameters() +// { +// // Set calibration parameters: +// // if calibration database exists, they are read from database, +// // otherwise, they are taken from digitizer. +// // +// // It is a user responsilibity to open CDB before reconstruction, +// // for example: +// // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); + +// //Check if calibration is stored in data base + +// if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet())) +// { +// AliCDBEntry *entry = (AliCDBEntry*) +// AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); +// if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); +// } + +// if(!fCalibData) +// printf("************* Calibration parameters not found in CDB! ****************"); +// // AliFatal("Calibration parameters not found in CDB!"); + + +// } + +// //____________________________________________________________________________ +// Float_t AliJetFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId) +// { + +// // Convert digitized amplitude into energy. +// // Calibration parameters are taken from calibration data base for raw data, +// // or from digitizer parameters for simulated data. + +// if(fCalibData){ + +// if (fGeom==0) +// printf("************* Did not get geometry from EMCALLoader ***************"); + +// Int_t iSupMod = -1; +// Int_t nModule = -1; +// Int_t nIphi = -1; +// Int_t nIeta = -1; +// Int_t iphi = -1; +// Int_t ieta = -1; + +// Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ; +// if(!bCell) { +// // fGeom->PrintGeometry(); +// Error("Calibrate()"," Wrong cell id number : %i", AbsId); +// assert(0); +// } + +// fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta); + +// fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi); +// fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi); + +// return -fADCpedestalECA + amp * fADCchannelECA ; + +// } +// else //Return energy with default parameters if calibration is not available +// return -fADCpedestalECA + amp * fADCchannelECA ; + +// } + + //_____________________________________________________________________________ Float_t AliJetFillUnitArrayEMCalDigits::EtaToTheta(Float_t arg) { diff --git a/JETAN/AliJetFillUnitArrayEMCalDigits.h b/JETAN/AliJetFillUnitArrayEMCalDigits.h index a127d4f7e39..a7abf204d5a 100644 --- a/JETAN/AliJetFillUnitArrayEMCalDigits.h +++ b/JETAN/AliJetFillUnitArrayEMCalDigits.h @@ -14,24 +14,22 @@ #include "TTask.h" #endif +class AliEMCALGeometry; class AliJetDummyGeo; class AliESDCaloCluster; -class AliEMCALCalibData; class AliJetReader; class AliJetESDReader; class TClonesArray; class TRefArray; class AliJetUnitArray; -//class AliESD; class AliESDEvent; class AliJetGrid; +class AliEMCALCalibData ; class AliJetFillUnitArrayEMCalDigits : public TTask { public: AliJetFillUnitArrayEMCalDigits(); - //PH AliJetFillUnitArrayEMCalDigits(Int_t event); - //PH AliJetFillUnitArrayEMCalDigits(AliESD *fESD); AliJetFillUnitArrayEMCalDigits(AliESDEvent *fESD); virtual ~AliJetFillUnitArrayEMCalDigits(); @@ -43,55 +41,66 @@ class AliJetFillUnitArrayEMCalDigits : public TTask void SetRefArray(TRefArray *refArray) {fRefArray = refArray;} void SetTPCGrid(AliJetGrid *grid) {fTPCGrid = grid;} void SetEMCalGrid(AliJetGrid *grid) {fEMCalGrid = grid;} - // void SetESD(AliESD *esd) {fESD = esd;} + void SetEleCorrection(Int_t flag = 1) {fECorrection = flag;} void SetESD(AliESDEvent *esd) {fESD = esd;} void SetInitMult(Int_t mult) {fNDigitEmcal = mult;} void SetInitMultCut(Int_t multcut) {fNDigitEmcalCut = multcut;} + void SetProcId(Bool_t id) {fProcId = id;} // Getter TClonesArray* GetUnitArray() {return fUnitArray;} TRefArray* GetRefArray() {return fRefArray;} Int_t GetMult() {return fNDigitEmcal;} Int_t GetMultCut() {return fNDigitEmcalCut;} + // For calibration => Under construction + // virtual Float_t Calibrate(Int_t amp, Int_t cellId) ; // Tranforms Amp to energy // Other void Exec(Option_t*); Float_t EtaToTheta(Float_t arg); protected: - AliESDEvent *fESD; // ESD - Int_t fNumUnits; // Number of units in the unit object array (same as num towers in EMCAL) - Float_t fEtaMinCal; // Define EMCAL acceptance in Eta - Float_t fEtaMaxCal; // Define EMCAL acceptance in Eta - Float_t fPhiMinCal; // Define EMCAL acceptance in Phi - Float_t fPhiMaxCal; // Define EMCAL acceptance in Phi - Int_t fNIn; // Number of Array filled in UnitArray - Int_t fOpt; // Detector to be used for jet reconstruction - Int_t fDebug; // Debug option - Int_t fNCEMCAL; // Number of clusters in EMCAL - Int_t fNCPHOS; // Number of clusters in PHOS - Int_t fNCCalo; // Number of cluster in EMCAL + PHOS calorimeters + AliESDEvent *fESD; // ESD + Int_t fNumUnits; // Number of units in the unit object array (same as num towers in EMCAL) + Float_t fEtaMinCal; // Define EMCAL acceptance in Eta + Float_t fEtaMaxCal; // Define EMCAL acceptance in Eta + Float_t fPhiMinCal; // Define EMCAL acceptance in Phi + Float_t fPhiMaxCal; // Define EMCAL acceptance in Phi + Int_t fNIn; // Number of Array filled in UnitArray + Int_t fOpt; // Detector to be used for jet reconstruction + Int_t fCluster; // Use all cells or cells in clusters for jet finding + Int_t fDebug; // Debug option + Int_t fNCEMCAL; // Number of clusters in EMCAL + Int_t fNCPHOS; // Number of clusters in PHOS + Int_t fNCCalo; // Number of cluster in EMCAL + PHOS calorimeters - AliJetGrid *fTPCGrid; // Define filled grid - AliJetGrid *fEMCalGrid; // Define filled grid + AliJetGrid *fTPCGrid; // Define filled grid + AliJetGrid *fEMCalGrid; // Define filled grid + Int_t fECorrection; // Electron correction flag - AliJetReaderHeader *fReaderHeader; // ReaderHeader - TClonesArray *fMomentumArray; // MomentumArray - TClonesArray *fUnitArray; // UnitArray - TRefArray *fRefArray; // UnitArray - AliJetDummyGeo *fGeom; // Set EMCal geometry + AliJetReaderHeader *fReaderHeader; // ReaderHeader + TClonesArray *fMomentumArray; // MomentumArray + TClonesArray *fUnitArray; // UnitArray + TRefArray *fRefArray; // UnitArray + Bool_t fProcId; // Bool_t for TProcessID synchronization + AliJetDummyGeo *fGeom; // Set EMCal geometry - AliESDCaloCluster *fClus; //! - Int_t fNDigitEmcal; //! - Int_t fNDigitEmcalCut; //! + AliESDCaloCluster *fClus; //! + Int_t fNDigitEmcal; //! + Int_t fNDigitEmcalCut; //! + //Calibration parameters... to be replaced by database + AliEMCALCalibData *fCalibData; //! Calibration database if aval + Float_t fADCchannelECA; // width of one ADC channel for EC section (GeV) + Float_t fADCpedestalECA; // pedestal of ADC for EC section (GeV) - private: AliJetFillUnitArrayEMCalDigits(const AliJetFillUnitArrayEMCalDigits &det); AliJetFillUnitArrayEMCalDigits &operator=(const AliJetFillUnitArrayEMCalDigits &det); void InitParameters(); - + // Under construction + // void GetCalibrationParameters(void) ; + ClassDef(AliJetFillUnitArrayEMCalDigits,1) // Fill Unit Array with tpc and/or emcal information }; diff --git a/JETAN/AliJetFillUnitArrayTracks.cxx b/JETAN/AliJetFillUnitArrayTracks.cxx index 18ab76657fd..b96ad08d45d 100644 --- a/JETAN/AliJetFillUnitArrayTracks.cxx +++ b/JETAN/AliJetFillUnitArrayTracks.cxx @@ -14,12 +14,13 @@ **************************************************************************/ -//---------------------------------------------------------------------- +//====================================================================== +// ***July 2006 // Fill Unit Array class -// Class used by AliJetESDReader to fill a UnitArray from the information -// extracted from the particle tracks +// Class used by AliJetESDReader to fill a UnitArray from the information extracted +// from the particle tracks // Author: magali.estienne@ires.in2p3.fr -//---------------------------------------------------------------------- +//====================================================================== // --- Standard library --- @@ -36,6 +37,7 @@ #include #include #include +#include // --- AliRoot header files --- #include "AliJetFinder.h" @@ -43,7 +45,6 @@ #include "AliJetReader.h" #include "AliJetESDReader.h" #include "AliJetESDReaderHeader.h" -//#include "AliESD.h" #include "AliESDEvent.h" #include "AliJetDummyGeo.h" #include "AliJetUnitArray.h" @@ -72,6 +73,7 @@ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks() fMomentumArray(0x0), fUnitArray(0x0), fRefArray(0x0), + fProcId(kFALSE), fTPCGrid(0x0), fEMCalGrid(0x0), fGeom(0x0), @@ -104,7 +106,7 @@ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks() } //_____________________________________________________________________________ -AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/) +AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* esd) : TTask("AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"), fNumUnits(0), fEtaMinCal(0), @@ -122,15 +124,16 @@ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/) fMomentumArray(0x0), fUnitArray(0x0), fRefArray(0x0), + fProcId(kFALSE), fTPCGrid(0x0), fEMCalGrid(0x0), fGeom(0x0), - fESD(0x0), + fESD(esd), fGrid0(0x0), fGrid1(0x0), fGrid2(0x0), fGrid3(0x0), - fGrid4(0x0), + fGrid4(0x0),// first: loop over tracks in ESD fNphi(0), fNeta(0), fPhi2(0x0), @@ -150,21 +153,124 @@ AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(AliESDEvent* /*esd*/) fPhiBinInEMCalAcc(0), fNbinPhi(0) { - // constructor + // constructor 2 +} + +//_____________________________________________________________________________ +AliJetFillUnitArrayTracks::AliJetFillUnitArrayTracks(const AliJetFillUnitArrayTracks &det) + : TTask(det),//"AliJetFillUnitArrayTracks","Fill Unit Array with tpc/its and emcal information"), + fNumUnits(det.fNumUnits), + fEtaMinCal(det.fEtaMinCal), + fEtaMaxCal(det.fEtaMaxCal), + fPhiMinCal(det.fPhiMinCal), + fPhiMaxCal(det.fPhiMaxCal), + fHadCorr(det.fHadCorr), + fHCorrection(det.fHCorrection), + fNTracks(det.fNTracks), + fNTracksCut(det.fNTracksCut), + fOpt(det.fOpt), + fDZ(det.fDZ), + fDebug(det.fDebug), + fReaderHeader(det.fReaderHeader), + fMomentumArray(det.fMomentumArray), + fUnitArray(det.fUnitArray), + fRefArray(det.fRefArray), + fProcId(det.fProcId), + fTPCGrid(det.fTPCGrid), + fEMCalGrid(det.fEMCalGrid), + fGeom(det.fGeom), + fESD(det.fESD), + fGrid0(det.fGrid0), + fGrid1(det.fGrid1), + fGrid2(det.fGrid2), + fGrid3(det.fGrid3), + fGrid4(det.fGrid4), + fNphi(det.fNphi), + fNeta(det.fNeta), + fPhi2(det.fPhi2), + fEta2(det.fEta2), + fPhi(det.fPhi), + fEta(det.fEta), + fIndex(det.fIndex), + fParams(det.fParams), + fGrid(det.fGrid), + fPhiMin(det.fPhiMin), + fPhiMax(det.fPhiMax), + fEtaMin(det.fEtaMin), + fEtaMax(det.fEtaMax), + fEtaBinInTPCAcc(det.fEtaBinInTPCAcc), + fPhiBinInTPCAcc(det.fPhiBinInTPCAcc), + fEtaBinInEMCalAcc(det.fEtaBinInEMCalAcc), + fPhiBinInEMCalAcc(det.fPhiBinInEMCalAcc), + fNbinPhi(det.fNbinPhi) +{ + // Copy constructor +} + +//------------------------------------------------------------------ +AliJetFillUnitArrayTracks& AliJetFillUnitArrayTracks::operator=(const AliJetFillUnitArrayTracks& other) +{ + // Assignment + + fNumUnits = other.fNumUnits; + fEtaMinCal = other.fEtaMinCal; + fEtaMaxCal = other.fEtaMaxCal; + fPhiMinCal = other.fPhiMinCal; + fPhiMaxCal = other.fPhiMaxCal; + fHadCorr = other.fHadCorr; + fHCorrection = other.fHCorrection; + fNTracks = other.fNTracks; + fNTracksCut = other.fNTracksCut; + fOpt = other.fOpt; + fDZ = other.fDZ; + fDebug = other.fDebug; + fReaderHeader = other.fReaderHeader; + fMomentumArray = other.fMomentumArray; + fUnitArray = other.fUnitArray; + fRefArray = other.fRefArray; + fProcId = other.fProcId; + fTPCGrid = other.fTPCGrid; + fEMCalGrid = other.fEMCalGrid; + fGeom = other.fGeom; + fESD = other.fESD; + fGrid0 = other.fGrid0; + fGrid1 = other.fGrid1; + fGrid2 = other.fGrid2; + fGrid3 = other.fGrid3; + fGrid4 = other.fGrid4; + fNphi = other.fNphi; + fNeta = other.fNeta; + fPhi2 = other.fPhi2; + fEta2 = other.fEta2; + fPhi = other.fPhi; + fEta = other.fEta; + fIndex = other.fIndex; + fParams = other.fParams; + fGrid = other.fGrid; + fPhiMin = other.fPhiMin; + fPhiMax = other.fPhiMax; + fEtaMin = other.fEtaMin; + fEtaMax = other.fEtaMax; + fEtaBinInTPCAcc = other.fEtaBinInTPCAcc; + fPhiBinInTPCAcc = other.fPhiBinInTPCAcc; + fEtaBinInEMCalAcc = other.fEtaBinInEMCalAcc; + fPhiBinInEMCalAcc = other.fPhiBinInEMCalAcc; + fNbinPhi = other.fNbinPhi; + + return (*this); } //____________________________________________________________________________ void AliJetFillUnitArrayTracks::InitParameters() { - // fHCorrection = 0; // For hadron correction - fHadCorr = 0; // For hadron correction - fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL - fDebug = fReaderHeader->GetDebug(); + fHadCorr = 0; // For hadron correction + fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL + //fDebug = fReaderHeader->GetDebug(); fEtaMinCal = fGeom->GetArm1EtaMin(); fEtaMaxCal = fGeom->GetArm1EtaMax(); fPhiMinCal = fGeom->GetArm1PhiMin(); - fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; // A verifier quelle doit etre la derniere valeur ! + fPhiMaxCal = fGeom->GetArm1PhiMax()-10.; fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin, fPhiMax,fEtaMin,fEtaMax); @@ -200,19 +306,24 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/) // fDebug = fReaderHeader->GetDebug(); - + // Set parameters InitParameters(); // get number of tracks in event (for the loop) Int_t goodTrack = 0; Int_t nt = 0; - Float_t pt, eta,phi; + Int_t nmax = 0; + Float_t pt,eta,phi; + // Int_t sflag = 0; TVector3 p3; + + if(fDebug>1) cout << "fESD in Fill array : " << fESD << endl; + nt = fESD->GetNumberOfTracks(); if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl; - - // temporary storage of signal and pt cut flag + + // temporary storage of signal and pt cut flag Int_t* sflag = new Int_t[nt]; Int_t* cflag = new Int_t[nt]; @@ -220,6 +331,8 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/) Float_t ptMin = fReaderHeader->GetPtCut(); Float_t etaMin = fReaderHeader->GetFiducialEtaMin(); Float_t etaMax = fReaderHeader->GetFiducialEtaMax(); + Float_t phiMin = fReaderHeader->GetFiducialPhiMin(); + Float_t phiMax = fReaderHeader->GetFiducialPhiMax(); fOpt = fReaderHeader->GetDetector(); fDZ = fReaderHeader->GetDZ(); @@ -234,52 +347,71 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/) fGrid = fTPCGrid->GetGridType(); - - //loop over tracks - for (Int_t it = 0; it < nt; it++) { - AliESDtrack *track = fESD->GetTrack(it); + // Loop over tracks + nmax = nt; + for (Int_t it = 0; it < nmax; it++) { + AliESDtrack *track; + track = fESD->GetTrack(it); UInt_t status = track->GetStatus(); Double_t mom[3]; track->GetPxPyPz(mom); + p3.SetXYZ(mom[0],mom[1],mom[2]); + vector vtmp(3); + vtmp[0] = mom[0]; vtmp[1] = mom[1]; vtmp[2] = mom[2]; pt = p3.Pt(); + Float_t mass = 0.; mass = track->GetMass(); - - if (((status & AliESDtrack::kITSrefit) == 0) || - ((status & AliESDtrack::kTPCrefit) == 0)) continue; // quality check + + if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check + if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly() && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check eta = p3.Eta(); phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi()); - + if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut + if ( (phi > phiMax) || (phi < phiMin)) continue; // checking phi cut + + sflag[goodTrack]=0; + if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; + cflag[goodTrack]=0; + if (pt > ptMin) cflag[goodTrack]=1; // pt cut + + if(fDebug>20) cout << "In Pythia: Track:" << it << ", eta: " << eta << ", phi: " << phi << ", pt: " << pt << endl; - // sflag -> not yet implemented !!!! - if(fGrid==0) { // Only TPC filled from its grid in its total acceptance Int_t idTPC = fTPCGrid->GetIndex(phi,eta); Bool_t ok = kFALSE; - + Bool_t ref = kFALSE; + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1); + uArray->SetUnitTrackID(it); - + Float_t unitEnergy = 0.; unitEnergy = uArray->GetUnitEnergy(); + // nTracksTpcOnly is necessary to count the number of candidate cells + // but it doesn't give the real multiplicity -> it will be extracted + // in the jet finder and for that, it is necessary to fill the Px,Py,Pz + // information for each tracks stored in a given unitcell if(unitEnergy==0.){ nTracksTpcOnly++; ok = kTRUE; + ref = kTRUE; } + // Fill energy in TPC acceptance uArray->SetUnitEnergy(unitEnergy + pt); - uArray->SetUnitPxPyPz(mom); - uArray->SetUnitMass(mass); + uArray->SetUnitPxPyPz(kFALSE,vtmp); + uArray->SetUnitMass(mass); // Pt cut flag if(uArray->GetUnitEnergy()0) { - uArray->SetUnitDetectorFlag(kAll); - } - if(uArray->GetUnitEnergy()>0){ + // Signal flag + if(sflag[goodTrack] == 1) { + uArray->SetUnitSignalFlag(kGood); + uArray->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray->GetUnitEnergy()>0 && ref){ + if(!fProcId) { + fProcId = kTRUE; + // delete fRefArray; + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + } fRefArray->Add(uArray); } - - sflag[goodTrack]=0; - if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; - cflag[goodTrack]=0; - if (pt > ptMin) cflag[goodTrack]=1; // pt cut - goodTrack++; - + } if(fGrid==1) { Int_t nElements = fTPCGrid->GetNEntries(); - // Fill track information in EMCAL acceptance if((eta >= fEtaMin && eta <= fEtaMax) && (phi >= fPhiMin && phi <= fPhiMax))// && { - // Include dead-zones if(fDZ) { @@ -329,173 +460,249 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/) Int_t n1 = fGrid1->GetNEntries(); Int_t n2 = fGrid2->GetNEntries(); Int_t n3 = fGrid3->GetNEntries(); - + if(phi >= phimin0 && phi <= phimax0){ Int_t id0 = fGrid0->GetIndex(phi,eta)-1; AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements); + uArray0->SetUnitTrackID(it); + Float_t uEnergy0 = uArray0->GetUnitEnergy(); Bool_t ok0 = kFALSE; + Bool_t ref0 = kFALSE; if(uEnergy0==0.){ nTracksEmcalDZ++; ok0 = kTRUE; + ref0 = kTRUE; } uArray0->SetUnitEnergy(uEnergy0+pt); + uArray0->SetUnitPxPyPz(kFALSE,vtmp); if(uArray0->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); else { uArray0->SetUnitCutFlag(kPtHigher); if(ok0) nTracksEmcalDZCut++; } - if(uArray0->GetUnitEnergy()>0) + if(sflag[goodTrack] == 1) { + uArray0->SetUnitSignalFlag(kGood); + uArray0->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray0->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray0->GetUnitEnergy()>0 && ref0){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0)); + fProcId = kTRUE; + } fRefArray->Add(uArray0); + } } if(phi >= phimin1 && phi <= phimax1){ Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0; AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements); uArray1->SetUnitTrackID(it); + Float_t uEnergy1 = uArray1->GetUnitEnergy(); Bool_t ok1 = kFALSE; + Bool_t ref1 = kFALSE; if(uEnergy1==0.){ nTracksEmcalDZ++; ok1 = kTRUE; + ref1 = kTRUE; } uArray1->SetUnitEnergy(uEnergy1+pt); + uArray1->SetUnitPxPyPz(kFALSE,vtmp); if(uArray1->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); else { uArray1->SetUnitCutFlag(kPtHigher); if(ok1) nTracksEmcalDZCut++; } - if(uArray1->GetUnitEnergy()>0) + if(sflag[goodTrack] == 1) { + uArray1->SetUnitSignalFlag(kGood); + uArray1->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray1->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray1->GetUnitEnergy()>0 && ref1){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1)); + fProcId = kTRUE; + } fRefArray->Add(uArray1); + } } if(phi >= phimin2 && phi <= phimax2){ Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1; AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements); uArray2->SetUnitTrackID(it); + Float_t uEnergy2 = uArray2->GetUnitEnergy(); Bool_t ok2 = kFALSE; + Bool_t ref2 = kFALSE; if(uEnergy2==0.){ nTracksEmcalDZ++; ok2 = kTRUE; + ref2 = kTRUE; } uArray2->SetUnitEnergy(uEnergy2+pt); + uArray2->SetUnitPxPyPz(kFALSE,vtmp); if(uArray2->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); else { uArray2->SetUnitCutFlag(kPtHigher); if(ok2) nTracksEmcalDZCut++; } - if(uArray2->GetUnitEnergy()>0) + if(sflag[goodTrack] == 1) { + uArray2->SetUnitSignalFlag(kGood); + uArray2->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray2->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray2->GetUnitEnergy()>0 && ref2){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2)); + fProcId = kTRUE; + } fRefArray->Add(uArray2); + } } if(phi >= phimin3 && phi <= phimax3){ Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2; AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements); uArray3->SetUnitTrackID(it); + Float_t uEnergy3 = uArray3->GetUnitEnergy(); Bool_t ok3 = kFALSE; + Bool_t ref3 = kFALSE; if(uEnergy3==0.){ nTracksEmcalDZ++; ok3 = kTRUE; + ref3 = kTRUE; } uArray3->SetUnitEnergy(uEnergy3+pt); + uArray3->SetUnitPxPyPz(kFALSE,vtmp); if(uArray3->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); else { uArray3->SetUnitCutFlag(kPtHigher); if(ok3) nTracksEmcalDZCut++; } - if(uArray3->GetUnitEnergy()>0) + if(sflag[goodTrack] == 1) { + uArray3->SetUnitSignalFlag(kGood); + uArray3->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray3->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray3->GetUnitEnergy()>0 && ref3){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3)); + fProcId = kTRUE; + } fRefArray->Add(uArray3); + } } if(phi >= phimin4 && phi <= phimax4){ Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3; AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements); uArray4->SetUnitTrackID(it); + Float_t uEnergy4 = uArray4->GetUnitEnergy(); Bool_t ok4 = kFALSE; + Bool_t ref4 = kFALSE; if(uEnergy4==0.){ nTracksEmcalDZ++; ok4 = kTRUE; + ref4 = kTRUE; } uArray4->SetUnitEnergy(uEnergy4+pt); + uArray4->SetUnitPxPyPz(kFALSE,vtmp); if(uArray4->GetUnitEnergy()SetUnitCutFlag(kPtSmaller); else { uArray4->SetUnitCutFlag(kPtHigher); if(ok4) nTracksEmcalDZCut++; } - if(uArray4->GetUnitEnergy()>0) + if(sflag[goodTrack] == 1) { + uArray4->SetUnitSignalFlag(kGood); + uArray4->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray4->SetUnitSignalFlagC(kFALSE,kBad); + + if(uArray4->GetUnitEnergy()>0 && ref4){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4)); + fProcId = kTRUE; + } fRefArray->Add(uArray4); + } } } // end fDZ Int_t towerID = 0; fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID); - if(towerID==-1) continue; - + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID); uArray->SetUnitTrackID(it); + Float_t unitEnergy = uArray->GetUnitEnergy(); Bool_t ok = kFALSE; + Bool_t ref = kFALSE; if(unitEnergy==0.){ nTracksEmcal++; ok=kTRUE; + ref=kTRUE; } // Do Hadron Correction // Parametrization to be added if (fHCorrection != 0) { - // Float_t hCEnergy = fHadCorr->GetEnergy(enT[i], (Double_t)etaT[i]); Float_t hCEnergy = fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta); unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta))); } //end Hadron Correction loop - + uArray->SetUnitEnergy(unitEnergy + pt); - + uArray->SetUnitPxPyPz(kFALSE,vtmp); // Put a pt cut flag - if(uArray->GetUnitEnergy()GetUnitEnergy()SetUnitCutFlag(kPtSmaller); + } else { uArray->SetUnitCutFlag(kPtHigher); if(ok) nTracksEmcalCut++; } - // Detector flag - if(unitEnergy > 0) - uArray->SetUnitDetectorFlag(kAll); - if(uArray->GetUnitEnergy()>0) - fRefArray->Add(uArray); - - sflag[goodTrack]=0; - if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; - cflag[goodTrack]=0; - if (pt > ptMin) cflag[goodTrack]=1; // pt cut - goodTrack++; + // Signal flag + if(sflag[goodTrack] == 1) { + uArray->SetUnitSignalFlag(kGood); + uArray->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray->SetUnitSignalFlagC(kFALSE,kBad); + if(uArray->GetUnitEnergy()>0 && ref){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + } } // end loop on EMCal acceptance cut + tracks quality else{ // Outside EMCal acceptance - // Int_t idTPC = GetIndexFromEtaPhi(etaT[i],phiT[i]); Int_t idTPC = fTPCGrid->GetIndex(phi,eta); - + AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC); uArray->SetUnitTrackID(it); Float_t unitEnergy2 = uArray->GetUnitEnergy(); // check if fNumUnits or fNumUnits-1 - Bool_t ok2 = kFALSE; + Bool_t okout = kFALSE; + Bool_t refout = kFALSE; if(unitEnergy2==0.){ nTracksTpc++; - ok2=kTRUE; + okout=kTRUE; + refout=kTRUE; } // Fill energy outside emcal acceptance uArray->SetUnitEnergy(unitEnergy2 + pt); + uArray->SetUnitPxPyPz(kFALSE,vtmp); // Pt cut flag if(uArray->GetUnitEnergy()SetUnitCutFlag(kPtHigher); - if(ok2) nTracksTpcCut++; + if(okout) nTracksTpcCut++; } - // Detector flag - if(unitEnergy2 > 0) - uArray->SetUnitDetectorFlag(kTpc); - if(uArray->GetUnitEnergy()>0) - fRefArray->Add(uArray); - - sflag[goodTrack]=0; - if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1; - cflag[goodTrack]=0; - if (pt > ptMin) cflag[goodTrack]=1; // pt cut - goodTrack++; + // Signal flag + if(sflag[goodTrack] == 1) { + uArray->SetUnitSignalFlag(kGood); + uArray->SetUnitSignalFlagC(kFALSE,kGood); + } else uArray->SetUnitSignalFlagC(kFALSE,kBad); + if(uArray->GetUnitEnergy()>0 && refout){ + if(!fProcId){ + new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray)); + fProcId = kTRUE; + } + fRefArray->Add(uArray); + } } } // end fGrid==1 + + goodTrack++; + } // end loop on entries (tpc tracks) -// // set the signal flags -// fSignalFlag.Set(goodTrack,sflag); -// fCutFlag.Set(goodTrack,cflag); + if(fDebug>0) + { + cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl; + cout << "goodTracks: " << goodTrack << endl; + } -// delete sflag; -// delete cflag; + delete[] sflag; + delete[] cflag; - // } // end loop on entries (tpc tracks) - if(fGrid==0) { fNTracks = nTracksTpcOnly; fNTracksCut = nTracksTpcOnlyCut; @@ -546,54 +757,53 @@ void AliJetFillUnitArrayTracks::Exec(Option_t* /*option*/) cout << "fNTracksCut : " << fNTracksCut << endl; } } - + } //__________________________________________________________ void AliJetFillUnitArrayTracks::GetEtaPhiFromIndex(Int_t index, Float_t &eta, Float_t &phi) { -// -// Obtain the index from eta and phi -// - for(Int_t j=0; jAt(i); - phi = fPhi2->At(j); - } - } - - // TPC-EMCAL grid - //------------------------------------- - Int_t ii = 0; - if(i==0) ii = 0; - if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i; - if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && iAt(i); - phi = fPhi2->At(j); - } - - if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) && - ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) { - if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);} - else eta = fEta2->At(i); - phi = fPhi2->At(j); - } - - if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) { - eta = fEta2->At(i); - phi = fPhi2->At(j); + for(Int_t j=0; jAt(i); + phi = fPhi2->At(j); + } + } + + // TPC-EMCAL grid + //------------------------------------- + Int_t ii = 0; + if(i==0) ii = 0; + if(i>0 && i<(fEtaBinInTPCAcc-fEtaBinInEMCalAcc)/2) ii = i; + if(i>=(fEtaBinInTPCAcc+fEtaBinInEMCalAcc)/2 && iAt(i); + phi = fPhi2->At(j); + } + + if((j>=(fNbinPhi+1) && j<(fNbinPhi+1+fPhiBinInEMCalAcc)) && + ((fNbinPhi+1)*(fNeta+1) + (j-fNbinPhi-1)*(fEtaBinInTPCAcc-fEtaBinInEMCalAcc) + ii)== index ) { + if(ii==0) {Int_t ind = 0; eta = fEta2->At(ind);} + else eta = fEta2->At(i); + phi = fPhi2->At(j); + } + + if(j>=(fNbinPhi+1+fPhiBinInEMCalAcc) && ((fNbinPhi+1)*(fNeta+1)+fPhiBinInEMCalAcc*((fEtaBinInTPCAcc-fEtaBinInEMCalAcc))+(j-(fNbinPhi+1+fPhiBinInEMCalAcc))*(fNeta+1)+i == index)) { + eta = fEta2->At(i); + phi = fPhi2->At(j); + } + } } - } } - } } diff --git a/JETAN/AliJetFillUnitArrayTracks.h b/JETAN/AliJetFillUnitArrayTracks.h index 6e15f312b63..6f8db05018c 100644 --- a/JETAN/AliJetFillUnitArrayTracks.h +++ b/JETAN/AliJetFillUnitArrayTracks.h @@ -17,12 +17,15 @@ #include #include +class AliEMCALGeometry; +class AliJetDummyGeo; class AliJetHadronCorrection; class AliJetReader; class AliJetESDReader; class TClonesArray; class TRefArray; class AliJetGrid; +class AliEMCALGeometry; class AliJetDummyGeo; class AliESD; class AliESDEvent; @@ -31,49 +34,48 @@ class AliJetFillUnitArrayTracks : public TTask { public: AliJetFillUnitArrayTracks(); - //PH AliJetFillUnitArrayTracks(AliESD *fESD); AliJetFillUnitArrayTracks(AliESDEvent *fESD); virtual ~AliJetFillUnitArrayTracks(); // Setter - void SetReaderHeader(AliJetReaderHeader *readerHeader) {fReaderHeader = readerHeader;} - void SetGeom(AliJetDummyGeo *geom) {fGeom = geom;} - void SetMomentumArray(TClonesArray *momentumArray) {fMomentumArray = momentumArray;} - void SetUnitArray(TClonesArray *unitArray) {fUnitArray = unitArray;} - void SetRefArray(TRefArray *refArray) {fRefArray = refArray;} - void SetHadCorrection(Int_t flag = 1) {fHCorrection = flag;} - void SetHadCorrector(AliJetHadronCorrectionv1* corr) {fHadCorr = corr;} - void SetTPCGrid(AliJetGrid *grid) {fTPCGrid = grid;} - void SetEMCalGrid(AliJetGrid *grid) {fEMCalGrid = grid;} - //PH void SetGrid(Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax); - // void SetESD(AliESD *esd) {fESD = esd;} - void SetESD(AliESDEvent *esd) {fESD = esd;} - void SetGrid0(AliJetGrid *grid0){fGrid0 = grid0;} - void SetGrid1(AliJetGrid *grid1){fGrid1 = grid1;} - void SetGrid2(AliJetGrid *grid2){fGrid2 = grid2;} - void SetGrid3(AliJetGrid *grid3){fGrid3 = grid3;} - void SetGrid4(AliJetGrid *grid4){fGrid4 = grid4;} + void SetReaderHeader(AliJetReaderHeader *readerHeader) {fReaderHeader = readerHeader;} + void SetGeom(AliJetDummyGeo *geom) {fGeom = geom;} + void SetMomentumArray(TClonesArray *momentumArray) {fMomentumArray = momentumArray;} + void SetUnitArray(TClonesArray *unitArray) {fUnitArray = unitArray;} + void SetRefArray(TRefArray *refArray) {fRefArray = refArray;} + void SetHadCorrection(Int_t flag = 1) {fHCorrection = flag;} + void SetHadCorrector(AliJetHadronCorrectionv1* corr) {fHadCorr = corr;} + void SetTPCGrid(AliJetGrid *grid) {fTPCGrid = grid;} + void SetEMCalGrid(AliJetGrid *grid) {fEMCalGrid = grid;} + void SetGrid(Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax); + void SetESD(AliESDEvent *esd) {fESD = esd;} + void SetGrid0(AliJetGrid *grid0){fGrid0 = grid0;} + void SetGrid1(AliJetGrid *grid1){fGrid1 = grid1;} + void SetGrid2(AliJetGrid *grid2){fGrid2 = grid2;} + void SetGrid3(AliJetGrid *grid3){fGrid3 = grid3;} + void SetGrid4(AliJetGrid *grid4){fGrid4 = grid4;} + void SetProcId(Bool_t id) {fProcId = id;} + Bool_t GetProcId() {return fProcId;} // Getter TClonesArray* GetUnitArray() {return fUnitArray;} TRefArray* GetRefArray() {return fRefArray;} - // Int_t GetIndexFromEtaPhi(Double_t eta,Double_t phi) const; void GetEtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi); - Int_t GetNeta() const {return fNeta;} - Int_t GetNphi() const {return fNphi;} - Int_t GetHadCorrection() const {return fHCorrection;} - Int_t GetMult() const {return fNTracks;} - Int_t GetMultCut() const {return fNTracksCut;} + Int_t GetNeta() {return fNeta;} + Int_t GetNphi() {return fNphi;} + Int_t GetHadCorrection() {return fHCorrection;} + Int_t GetMult() {return fNTracks;} + Int_t GetMultCut() {return fNTracksCut;} void Exec(Option_t*); protected: - Int_t fNumUnits; // Number of units in the unit object array (same as num towers in EMCAL) - Float_t fEtaMinCal; // Define EMCal acceptance in Eta - Float_t fEtaMaxCal; // Define EMCal acceptance in Eta - Float_t fPhiMinCal; // Define EMCal acceptance in Phi - Float_t fPhiMaxCal; // Define EMCal acceptance in Phi + Int_t fNumUnits; // Number of units in the unit object array (same as num towers in EMCAL) + Float_t fEtaMinCal; // Define EMCal acceptance in Eta + Float_t fEtaMaxCal; // Define EMCal acceptance in Eta + Float_t fPhiMinCal; // Define EMCal acceptance in Phi + Float_t fPhiMaxCal; // Define EMCal acceptance in Phi AliJetHadronCorrectionv1 *fHadCorr; // Pointer to Hadron Correction Object - Int_t fHCorrection; // Hadron correction flag + Int_t fHCorrection; // Hadron correction flag Int_t fNTracks; // Number of tracks stored in UnitArray Int_t fNTracksCut; // Number of tracks stored in UnitArray with a pt cut Int_t fOpt; // Detector to be used for jet reconstruction @@ -84,40 +86,40 @@ class AliJetFillUnitArrayTracks : public TTask TClonesArray *fMomentumArray; // MomentumArray TClonesArray *fUnitArray; // UnitArray TRefArray *fRefArray; // UnitArray + Bool_t fProcId; // Bool_t for TProcessID synchronization AliJetGrid *fTPCGrid; // Define filled grid AliJetGrid *fEMCalGrid; // Define filled grid AliJetDummyGeo *fGeom; // Define EMCal geometry AliESDEvent *fESD; // ESD - AliJetGrid *fGrid0; // Pointer to Grid 1 - AliJetGrid *fGrid1; // Pointer to Grid 2 - AliJetGrid *fGrid2; // Pointer to Grid 3 - AliJetGrid *fGrid3; // Pointer to Grid 4 - AliJetGrid *fGrid4; // Pointer to Grid 5 + AliJetGrid *fGrid0; // Grid used for dead zones definition + AliJetGrid *fGrid1; // Grid used for dead zones definition + AliJetGrid *fGrid2; // Grid used for dead zones definition + AliJetGrid *fGrid3; // Grid used for dead zones definition + AliJetGrid *fGrid4; // Grid used for dead zones definition - Int_t fNphi; // number of points in the grid: phi - Int_t fNeta; // " eta - TArrayD* fPhi2; // grid points in phi - TArrayD* fEta2; // grid points in eta - TArrayD* fPhi; // grid points in phi - TArrayD* fEta; // grid points in eta - TMatrixD* fIndex; // grid points in (phi,eta) - TMatrixD* fParams; // matrix of parameters in the grid points - Int_t fGrid; // Select the grid acceptance you want to fill - // 0 = TPC acceptance, 1 = TPC-EMCal acceptance - Float_t fPhiMin; // minimum phi - Float_t fPhiMax; // maximum phi - Float_t fEtaMin; // minimum eta - Float_t fEtaMax; // maximum eta - Int_t fEtaBinInTPCAcc; // eta bins in tpc acceptance - Int_t fPhiBinInTPCAcc; // phi bins in tpc acceptance - Int_t fEtaBinInEMCalAcc; // eta bins in emcal acceptance - Int_t fPhiBinInEMCalAcc; // phi bins in emcal acceptance - Int_t fNbinPhi; // number of phi bins + Int_t fNphi; // number of points in the grid: phi + Int_t fNeta; // " eta + TArrayD *fPhi2; // grid points in phi + TArrayD *fEta2; // grid points in eta + TArrayD *fPhi; // grid points in phi + TArrayD *fEta; // grid points in eta + TMatrixD *fIndex; // grid points in (phi,eta) + TMatrixD *fParams; // matrix of parameters in the grid points + Int_t fGrid; // Select the grid acceptance you want to fill + // 0 = TPC acceptance, 1 = TPC-EMCal acceptance + Float_t fPhiMin; + Float_t fPhiMax; + Float_t fEtaMin; + Float_t fEtaMax; + Int_t fEtaBinInTPCAcc; + Int_t fPhiBinInTPCAcc; + Int_t fEtaBinInEMCalAcc; + Int_t fPhiBinInEMCalAcc; + Int_t fNbinPhi; private: AliJetFillUnitArrayTracks(const AliJetFillUnitArrayTracks &det); AliJetFillUnitArrayTracks &operator=(const AliJetFillUnitArrayTracks &det); - // void SetEMCALGeometry(); void InitParameters(); ClassDef(AliJetFillUnitArrayTracks,1) // Fill Unit Array with tpc and/or emcal information diff --git a/JETAN/AliJetFinder.cxx b/JETAN/AliJetFinder.cxx index d8cb7ba9e12..25239894afb 100644 --- a/JETAN/AliJetFinder.cxx +++ b/JETAN/AliJetFinder.cxx @@ -20,11 +20,13 @@ // manages the search for jets // Authors: jgcn@mda.cinvestav.mx // andreas.morsch@cern.ch +// magali.estienne@subatech.in2p3.fr //--------------------------------------------------------------------- #include #include #include +#include #include "AliJetFinder.h" #include "AliJet.h" @@ -34,6 +36,7 @@ #include "AliJetControlPlots.h" #include "AliLeading.h" #include "AliAODEvent.h" +#include "AliJetUnitArray.h" ClassImp(AliJetFinder) @@ -51,7 +54,10 @@ AliJetFinder::AliJetFinder(): fOut(0x0) { + // // Constructor + // + fJets = new AliJet(); fGenJets = new AliJet(); fLeading = new AliLeading(); @@ -59,11 +65,13 @@ AliJetFinder::AliJetFinder(): } //////////////////////////////////////////////////////////////////////// - AliJetFinder::~AliJetFinder() { - // destructor - // here reset and delete jets + // + // Destructor + // + + // Reset and delete jets fJets->ClearJets(); delete fJets; fGenJets->ClearJets(); @@ -76,13 +84,9 @@ AliJetFinder::~AliJetFinder() delete fOut; // reset and delete control plots if (fPlots) delete fPlots; - if ( fLeading ) - delete fLeading; - fLeading = NULL; } //////////////////////////////////////////////////////////////////////// - void AliJetFinder::SetOutputFile(const char */*name*/) { // opens output file @@ -90,7 +94,6 @@ void AliJetFinder::SetOutputFile(const char */*name*/) } //////////////////////////////////////////////////////////////////////// - void AliJetFinder::PrintJets() { // Print jet information @@ -101,7 +104,6 @@ void AliJetFinder::PrintJets() } //////////////////////////////////////////////////////////////////////// - void AliJetFinder::SetPlotMode(Bool_t b) { // Sets the plotting mode @@ -121,7 +123,6 @@ TTree* AliJetFinder::MakeTreeJ(char* name) } //////////////////////////////////////////////////////////////////////// - void AliJetFinder::WriteRHeaderToFile() { // write reader header @@ -130,8 +131,6 @@ void AliJetFinder::WriteRHeaderToFile() } //////////////////////////////////////////////////////////////////////// - - void AliJetFinder::Run() { // Do some initialization @@ -158,11 +157,11 @@ void AliJetFinder::Run() if (option == 0) { // TPC with fMomentumArray if(debug > 1) - printf("In FindJetsTPC() routine: find jets with fMomentumArray !!!\n"); - FindJetsTPC(); + printf("In FindJetsC() routine: find jets with fMomentumArray !!!\n"); + FindJetsC(); } else { - if(debug > 1) printf("In FindJets() routine: find jets with fUnitArray !!!\n"); - FindJets(); + if(debug > 1) printf("In FindJets() routine: find jets with fUnitArray !!!\n"); + FindJets(); } if (fOut) { fOut->cd(); @@ -191,12 +190,14 @@ void AliJetFinder::Run() // The following methods have been added to allow for event steering from the outside // +//////////////////////////////////////////////////////////////////////// void AliJetFinder::ConnectTree(TTree* tree, TObject* data) { // Connect the input file fReader->ConnectTree(tree, data); } +//////////////////////////////////////////////////////////////////////// void AliJetFinder::WriteHeaders() { // Write the Headers @@ -206,27 +207,84 @@ void AliJetFinder::WriteHeaders() f->Close(); } - +//////////////////////////////////////////////////////////////////////// Bool_t AliJetFinder::ProcessEvent() { -// -// Process one event -// - Bool_t ok = fReader->FillMomentumArray(); - if (!ok) return kFALSE; - - // Leading particles - fLeading->FindLeading(fReader); - // Jets - FindJets(); - - if (fPlots) fPlots->FillHistos(fJets); - fLeading->Reset(); - fGenJets->ClearJets(); - Reset(); - return kTRUE; + // + // Process one event + // Charged only jets + // + + Bool_t ok = fReader->FillMomentumArray(); + if (!ok) return kFALSE; + + // Leading particles + fLeading->FindLeading(fReader); + // Jets + FindJets(); // V1 + // FindJetsC(); // V2 + + if (fPlots) fPlots->FillHistos(fJets); + fLeading->Reset(); + fGenJets->ClearJets(); + Reset(); + return kTRUE; } +//////////////////////////////////////////////////////////////////////// +Bool_t AliJetFinder::ProcessEvent2() +{ + // + // Process one event + // Charged only or charged+neutral jets + // + + TRefArray* ref = new TRefArray(); + Bool_t procid = kFALSE; + Bool_t ok = fReader->ExecTasks(procid,ref); + + // Delete reference pointer + if (!ok) {delete ref; return kFALSE;} + + // Leading particles + fLeading->FindLeading(fReader); + // Jets + FindJets(); + + Int_t nEntRef = ref->GetEntries(); + vector vtmp; + + for(Int_t i=0; iAt(i))->SetUnitTrackID(0); + ((AliJetUnitArray*)ref->At(i))->SetUnitEnergy(0.); + ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag(kPtSmaller); + ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag2(kPtSmaller); + ((AliJetUnitArray*)ref->At(i))->SetUnitPxPyPz(kTRUE,vtmp); + ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlag(kBad); + ((AliJetUnitArray*)ref->At(i))->SetUnitDetectorFlag(kTpc); + ((AliJetUnitArray*)ref->At(i))->SetUnitFlag(kOutJet); + + // Reset process ID + AliJetUnitArray* uA = (AliJetUnitArray*)ref->At(i); + uA->ResetBit(kIsReferenced); + uA->SetUniqueID(0); + } + + // Delete the reference pointer + ref->Delete(); + delete ref; + + if (fPlots) fPlots->FillHistos(fJets); + fLeading->Reset(); + fGenJets->ClearJets(); + Reset(); + + return kTRUE; +} + +//////////////////////////////////////////////////////////////////////// void AliJetFinder::FinishRun() { // Finish a run @@ -244,21 +302,25 @@ void AliJetFinder::FinishRun() } } +//////////////////////////////////////////////////////////////////////// void AliJetFinder::AddJet(AliAODJet p) { - // Add new jet to the list +// Add new jet to the list new ((*fAODjets)[fNAODjets++]) AliAODJet(p); } +//////////////////////////////////////////////////////////////////////// void AliJetFinder::ConnectAOD(AliAODEvent* aod) { // Connect to the AOD fAODjets = aod->GetJets(); } +//////////////////////////////////////////////////////////////////////// void AliJetFinder::ConnectAODNonStd(AliAODEvent* aod,const char *bname) { fAODjets = dynamic_cast(aod->FindListObject(bname)); // how is this is reset? Cleared? } + diff --git a/JETAN/AliJetFinder.h b/JETAN/AliJetFinder.h index ec62bd131b9..a782d7c8e3b 100755 --- a/JETAN/AliJetFinder.h +++ b/JETAN/AliJetFinder.h @@ -9,6 +9,7 @@ // manages the search for jets // Authors: jgcn@mda.cinvestav.mx // andreas.morsch@cern.ch +// magali.estienne@subatech.in2p3.fr //--------------------------------------------------------------------- #include @@ -16,7 +17,7 @@ #include "AliJetReader.h" class TFile; -class TTree; +class TChain; class AliJet; class AliJetHeader; class AliJetControlPlots; @@ -33,35 +34,38 @@ class AliJetFinder : public TObject // getters virtual AliJetReader *GetReader() {return fReader;} virtual AliJetHeader *GetHeader() {return fHeader;} - virtual AliJet *GetJets() {return fJets;} - virtual Bool_t GetPlotMode() const {return fPlotMode;} - virtual TFile* GetOutputFile() {return fOut;} + virtual AliJet *GetJets() {return fJets;} + virtual Bool_t GetPlotMode() const {return fPlotMode;} + virtual TFile *GetOutputFile() {return fOut;} // setters - virtual void SetPlotMode(Bool_t b); - virtual void SetOutputFile(const char *name="jets.root"); - virtual void SetJetReader(AliJetReader* r) {fReader=r;} - virtual void SetJetHeader(AliJetHeader* h) {fHeader=h;} + virtual void SetPlotMode(Bool_t b); + virtual void SetOutputFile(const char *name="jets.root"); + virtual void SetJetReader(AliJetReader* r) {fReader=r;} + virtual void SetJetHeader(AliJetHeader* h) {fHeader=h;} // others - virtual void AddJet(AliAODJet jet); - virtual void PrintJets(); - virtual void Run(); - virtual void WriteRHeaderToFile(); + virtual void AddJet(AliAODJet jet); + virtual void PrintJets(); + virtual void Run(); + virtual void WriteRHeaderToFile(); // the following have to be implemented for each specific finder - virtual void Init() {} - virtual void Reset() {fNAODjets = 0;} - virtual void FindJets() {} - virtual void FindJetsTPC(){} - virtual void WriteJHeaderToFile() { } + virtual void Init() {} + virtual void InitTask(TChain* /*tree*/) {} + virtual void Reset() {fNAODjets = 0;} + virtual void FindJets() {} + virtual void FindJetsC(){} + virtual void WriteJHeaderToFile() { } // some methods to allow steering from the outside - virtual Bool_t ProcessEvent(); - virtual void FinishRun(); - virtual void ConnectTree(TTree* tree, TObject* data); - virtual void ConnectAOD(AliAODEvent* aod); - virtual void ConnectAODNonStd(AliAODEvent* aod,const char* bname); - virtual TTree* MakeTreeJ(char* name); - virtual void WriteHeaders(); - virtual void WriteJetsToFile() {;} - virtual void TestJet() {;} + virtual Bool_t ProcessEvent(); + virtual Bool_t ProcessEvent2(); + virtual void FinishRun(); + virtual void ConnectTree(TTree* tree, TObject* data); + virtual void ConnectAOD(AliAODEvent* aod); + virtual void ConnectAODNonStd(AliAODEvent* aod,const char* bname); + virtual TTree *MakeTreeJ(char* name); + virtual void WriteHeaders(); + virtual void WriteJetsToFile() {;} + virtual void TestJet() {;} + protected: AliJetFinder(const AliJetFinder& rJetFinder); AliJetFinder& operator = (const AliJetFinder& rhsf); @@ -76,6 +80,7 @@ class AliJetFinder : public TObject Int_t fNAODjets; //! number of reconstructed jets AliJetControlPlots* fPlots; //! pointer to control plots TFile* fOut; //! output file + ClassDef(AliJetFinder,2) }; diff --git a/JETAN/AliJetGrid.cxx b/JETAN/AliJetGrid.cxx index 55b30229eb3..9f8e497aa1c 100644 --- a/JETAN/AliJetGrid.cxx +++ b/JETAN/AliJetGrid.cxx @@ -30,7 +30,7 @@ // > grid->SetMatrixIndexes(); // > grid->SetIndexIJ(); // -// Author : magali.estienne@ires.in2p3.fr +// Author : magali.estienne@subatech.in2p3.fr //========================================================================= // Standard headers @@ -45,70 +45,75 @@ ClassImp(AliJetGrid) - //__________________________________________________________ AliJetGrid::AliJetGrid(): - fGrid(0), - fNphi(0), - fNeta(0), - fPhi(0), - fEta(0), - fIndex(0), - fIndexI(0), - fIndexJ(0), - fPhiMin(0), - fPhiMax(0), - fEtaMin(0), - fEtaMax(0), - fEtaBinInTPCAcc(0), - fPhiBinInTPCAcc(0), - fEtaBinInEMCalAcc(0), - fPhiBinInEMCalAcc(0), - fNbinEta(0), - fNbinPhi(0), - fMaxPhi(0), - fMinPhi(0), - fMaxEta(0), - fMinEta(0), - fDebug(1) - { + fGrid(0), + fNphi(0), + fNeta(0), + fPhi(0), + fEta(0), + fIndex(0), + fIndexI(0), + fIndexJ(0), + fPhiMin(0), + fPhiMax(0), + fEtaMin(0), + fEtaMax(0), + fEtaBinInTPCAcc(0), + fPhiBinInTPCAcc(0), + fEtaBinInEMCalAcc(0), + fPhiBinInEMCalAcc(0), + fNbinEta(0), + fNbinPhi(0), + fMaxPhi(0), + fMinPhi(0), + fMaxEta(0), + fMinEta(0), + fDebug(1) +{ // Default constructor } //__________________________________________________________ AliJetGrid::AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t phiMax,Double_t etaMin,Double_t etaMax): - fGrid(0), - fNphi(nphi), - fNeta(neta), - fPhi(0), - fEta(0), - fIndex(0), - fIndexI(0), - fIndexJ(0), - fPhiMin(0), - fPhiMax(0), - fEtaMin(0), - fEtaMax(0), - fEtaBinInTPCAcc(0), - fPhiBinInTPCAcc(0), - fEtaBinInEMCalAcc(0), - fPhiBinInEMCalAcc(0), - fNbinEta(0), - fNbinPhi(0), - fMaxPhi(phiMax), - fMinPhi(phiMin), - fMaxEta(etaMax), - fMinEta(etaMin), - fDebug(1) - { + fGrid(0), + fNphi(nphi), + fNeta(neta), + fPhi(0), + fEta(0), + fIndex(0), + fIndexI(0), + fIndexJ(0), + fPhiMin(0), + fPhiMax(0), + fEtaMin(0), + fEtaMax(0), + fEtaBinInTPCAcc(0), + fPhiBinInTPCAcc(0), + fEtaBinInEMCalAcc(0), + fPhiBinInEMCalAcc(0), + fNbinEta(0), + fNbinPhi(0), + fMaxPhi(phiMax), + fMinPhi(phiMin), + fMaxEta(etaMax), + fMinEta(etaMin), + fDebug(1) +{ // Standard constructor fPhi = new TArrayD(fNphi+1); fEta = new TArrayD(fNeta+1); fIndexI = new TArrayI((fNeta+1)*(fNphi+1)+1); fIndexJ = new TArrayI((fNeta+1)*(fNphi+1)+1); - for(Int_t i=0; i 3){ for(Int_t i=0; i<(fNphi+1); i++) cout << (*fPhi)[i] << endl; @@ -120,30 +125,31 @@ AliJetGrid::AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t phiMax,Dou } //__________________________________________________________ -AliJetGrid::AliJetGrid(const AliJetGrid& grid) : TNamed(grid), - fGrid(grid.fGrid), - fNphi(grid.fNphi), - fNeta(grid.fNeta), - fPhi(0), - fEta(0), - fIndex(0), - fIndexI(grid.fIndexI), - fIndexJ(grid.fIndexJ), - fPhiMin(grid.fPhiMin), - fPhiMax(grid.fPhiMax), - fEtaMin(grid.fEtaMin), - fEtaMax(grid.fEtaMax), - fEtaBinInTPCAcc(grid.fEtaBinInTPCAcc), - fPhiBinInTPCAcc(grid.fPhiBinInTPCAcc), - fEtaBinInEMCalAcc(grid.fEtaBinInEMCalAcc), - fPhiBinInEMCalAcc(grid.fPhiBinInEMCalAcc), - fNbinEta(grid.fNbinEta), - fNbinPhi(grid.fNbinPhi), - fMaxPhi(grid.fMaxPhi), - fMinPhi(grid.fMinPhi), - fMaxEta(grid.fMaxEta), - fMinEta(grid.fMinEta), - fDebug(grid.fDebug) +AliJetGrid::AliJetGrid(const AliJetGrid& grid) : + TNamed(grid), + fGrid(grid.fGrid), + fNphi(grid.fNphi), + fNeta(grid.fNeta), + fPhi(0), + fEta(0), + fIndex(0), + fIndexI(grid.fIndexI), + fIndexJ(grid.fIndexJ), + fPhiMin(grid.fPhiMin), + fPhiMax(grid.fPhiMax), + fEtaMin(grid.fEtaMin), + fEtaMax(grid.fEtaMax), + fEtaBinInTPCAcc(grid.fEtaBinInTPCAcc), + fPhiBinInTPCAcc(grid.fPhiBinInTPCAcc), + fEtaBinInEMCalAcc(grid.fEtaBinInEMCalAcc), + fPhiBinInEMCalAcc(grid.fPhiBinInEMCalAcc), + fNbinEta(grid.fNbinEta), + fNbinPhi(grid.fNbinPhi), + fMaxPhi(grid.fMaxPhi), + fMinPhi(grid.fMinPhi), + fMaxEta(grid.fMaxEta), + fMinEta(grid.fMinEta), + fDebug(grid.fDebug) { // Copy constructor @@ -162,40 +168,40 @@ AliJetGrid::AliJetGrid(const AliJetGrid& grid) : TNamed(grid), AliJetGrid& AliJetGrid::operator=(const AliJetGrid& other) { -// Assignment - fGrid = other.fGrid; - fNphi = other.fNphi; - fNeta = other.fNeta; - fPhi = 0; - fEta = 0; - fIndex = 0; - fIndexI = other.fIndexI; - fIndexJ = other.fIndexJ; - fPhiMin = other.fPhiMin; - fPhiMax = other.fPhiMax; - fEtaMin = other.fEtaMin; - fEtaMax = other.fEtaMax; - fEtaBinInTPCAcc = other.fEtaBinInTPCAcc; - fPhiBinInTPCAcc = other.fPhiBinInTPCAcc; - fEtaBinInEMCalAcc = other.fEtaBinInEMCalAcc; - fPhiBinInEMCalAcc = other.fPhiBinInEMCalAcc; - fNbinEta = other.fNbinEta; - fNbinPhi = other.fNbinPhi; - fMaxPhi = other.fMaxPhi; - fMinPhi = other.fMinPhi; - fMaxEta = other.fMaxEta; - fMinEta = other.fMinEta; - fDebug = other.fDebug; - fPhi = new TArrayD(fNphi+1); - for(Int_t i=0; iAt(i); - fEta = new TArrayD(fNeta+1); - for(Int_t i=0; iAt(i); - - fIndex = new TMatrixD(fNphi+1,fNeta+1); - for(Int_t i=0; iAt(i); + fEta = new TArrayD(fNeta+1); + for(Int_t i=0; iAt(i); + + fIndex = new TMatrixD(fNphi+1,fNeta+1); + for(Int_t i=0; i @@ -22,8 +24,6 @@ class AliJetGrid : public TNamed { AliJetGrid(); AliJetGrid(Int_t nphi,Int_t neta,Double_t phiMin,Double_t etaMin,Double_t phiMax,Double_t etaMax); - AliJetGrid(const AliJetGrid& grid); - AliJetGrid& operator=(const AliJetGrid& other); virtual ~AliJetGrid(); // Getter @@ -51,12 +51,10 @@ class AliJetGrid : public TNamed { void GetEtaPhiFromIndex2(Int_t index, Float_t &phi, Float_t &eta); Int_t GetNEntries(); Int_t GetNEntries2(); - Int_t GetDeta() const {return static_cast((fEtaMax-fEtaMin)/fNeta); - if(fDebug>21) cout << "static_cast((fEtaMax-fEtaMin)/fNeta) : " << - static_cast((fEtaMax-fEtaMin)/fNeta);} - Int_t GetDphi() const {return static_cast((fPhiMax-fPhiMin)/fNphi); - if(fDebug>21) cout << "static_cast((fPhiMax-fPhiMin)/fNphi) : " << - static_cast((fPhiMax-fPhiMin)/fNphi);} + Int_t GetDeta() const {if(fNeta!=0) return static_cast((fEtaMax-fEtaMin)/fNeta); + else return static_cast(fEtaMax-fEtaMin);} + Int_t GetDphi() const {if(fNphi!=0) return static_cast((fPhiMax-fPhiMin)/fNphi); + else return static_cast(fPhiMax-fPhiMin);} Int_t GetGridType() const {return fGrid;} // Setter @@ -97,7 +95,12 @@ class AliJetGrid : public TNamed { Double_t fMinEta; // minimum eta Int_t fDebug; // debug flag - ClassDef(AliJetGrid,1) // Parameters used by AliTPCtrackerParam + protected: + AliJetGrid(const AliJetGrid& grid); + AliJetGrid& operator=(const AliJetGrid& other); + + + ClassDef(AliJetGrid,1) }; diff --git a/JETAN/AliJetReader.cxx b/JETAN/AliJetReader.cxx index 18e2005f5da..6c761199d44 100755 --- a/JETAN/AliJetReader.cxx +++ b/JETAN/AliJetReader.cxx @@ -37,24 +37,30 @@ ClassImp(AliJetReader) //////////////////////////////////////////////////////////////////////// + AliJetReader::AliJetReader(): // Constructor fChain(0), - fMomentumArray(new TClonesArray("TLorentzVector",2000)), + fTree(0), + fMomentumArray(new TClonesArray("TLorentzVector",4000)), fArrayMC(0), fFillUnitArray(new TTask("fillUnitArray","Fill unit array jet finder")), fESD(0), fReaderHeader(0), + fAliHeader(0), fSignalFlag(0), fCutFlag(0), fUnitArray(new TClonesArray("AliJetUnitArray",60000)), - fRefArray(new TRefArray()), fUnitArrayNoCuts(new TClonesArray("AliJetUnitArray",60000)), fArrayInitialised(0), fFillUAFromTracks(new AliJetFillUnitArrayTracks()), fFillUAFromEMCalDigits(new AliJetFillUnitArrayEMCalDigits()), fNumCandidate(0), - fNumCandidateCut(0) + fNumCandidateCut(0), + fHadronCorrector(0), + fHCorrection(0), + fECorrection(0), + fEFlag(kFALSE) { // Default constructor fSignalFlag = TArrayI(); @@ -66,36 +72,33 @@ AliJetReader::AliJetReader(): AliJetReader::~AliJetReader() { // Destructor - if (fMomentumArray) { - fMomentumArray->Delete(); - delete fMomentumArray; - } - - if (fUnitArray) { - fUnitArray->Delete(); - delete fUnitArray; - } - - if (fUnitArrayNoCuts) { - fUnitArrayNoCuts->Delete(); - delete fUnitArrayNoCuts; - } - - if (fFillUnitArray) { - delete fFillUnitArray; - } - - if (fArrayMC) { - fArrayMC->Delete(); - delete fArrayMC; - } + if (fMomentumArray) { + fMomentumArray->Delete(); + delete fMomentumArray; + } + + if (fUnitArray) { + fUnitArray->Delete(); + delete fUnitArray; + } + + if (fUnitArrayNoCuts) { + fUnitArrayNoCuts->Delete(); + delete fUnitArrayNoCuts; + } + + if (fFillUnitArray) { + fFillUnitArray->Delete(); + delete fFillUnitArray; + } + delete fArrayMC; + } //////////////////////////////////////////////////////////////////////// void AliJetReader::ClearArray() - { if (fMomentumArray) fMomentumArray->Clear(); if (fFillUnitArray) fFillUnitArray->Clear(); diff --git a/JETAN/AliJetReader.h b/JETAN/AliJetReader.h index a27a683f1d8..457f3886a29 100755 --- a/JETAN/AliJetReader.h +++ b/JETAN/AliJetReader.h @@ -7,15 +7,15 @@ // Jet reader base class // manages the reading of input for jet algorithms // Authors: jgcn@mda.cinvestav.mx -// Magali Estienne +// Magali Estienne #include #include #include class TTree; -class TTask; class TChain; +class TTask; class TClonesArray; class TRefArray; class AliJetReaderHeader; @@ -34,13 +34,13 @@ class AliJetReader : public TObject virtual ~AliJetReader(); // Getters - virtual TClonesArray* GetMomentumArray() const {return fMomentumArray;} - virtual TRefArray* GetReferences() const {return 0;} - virtual TClonesArray *GetUnitArray() const {return fUnitArray;} - virtual TRefArray *GetRefArray() const {return fRefArray;} + virtual TClonesArray* GetMomentumArray() const {return fMomentumArray;} + virtual TRefArray* GetReferences() const {return 0;} + virtual TClonesArray *GetUnitArray() const {return fUnitArray;} virtual TClonesArray *GetUnitArrayNoCuts() const {return fUnitArrayNoCuts;} virtual AliJetReaderHeader* GetReaderHeader() const {return fReaderHeader;} + virtual AliHeader *GetAliHeader() const {return fAliHeader;} virtual Int_t GetSignalFlag(Int_t i) const {return fSignalFlag[i];} virtual Int_t GetCutFlag(Int_t i) const {return fCutFlag[i];} virtual Int_t GetArrayInitialised() const {return fArrayInitialised;} @@ -49,24 +49,21 @@ class AliJetReader : public TObject // Setters virtual Bool_t FillMomentumArray() {return kTRUE;} + virtual Bool_t ReadEventLoader(Int_t) {return kTRUE;} virtual void FillUnitArrayFromTPCTracks(Int_t) {} // temporarily not used virtual void FillUnitArrayFromEMCALHits() {} // temporarily not used virtual void FillUnitArrayFromEMCALDigits(Int_t) {} // temporarily not used - virtual void FillUnitArrayFromEMCALClusters(Int_t) {} // temporarily not used virtual void InitUnitArray() {} virtual void InitParameters() {} - virtual void CreateTasks() {} - // virtual void ExecTasks(Int_t) {} - virtual Bool_t ExecTasks(Int_t) {return kTRUE;} - /* // Correction of hadronic energy - virtual void SetHadronCorrector(AliEMCALHadronCorrectionv1* corr) {fHadronCorrector = corr;} - virtual void SetHadronCorrection(Int_t flag = 1) {fHCorrection = flag;} */ + virtual void CreateTasks(TChain* /*tree*/) {} + virtual Bool_t ExecTasks(Bool_t /*procid*/, TRefArray* /*refArray*/) {return kFALSE;} + // Correction of hadronic energy + virtual void SetHadronCorrector(AliJetHadronCorrectionv1*) {;} + virtual void SetHadronCorrection(Int_t) {;} + virtual void SetElectronCorrection(Int_t) {;} virtual void SetReaderHeader(AliJetReaderHeader* header) - {fReaderHeader = header;} + {fReaderHeader = header;} virtual void SetESD(AliESDEvent* esd) { fESD = esd;} - // virtual Int_t SetNumCandidate(Int_t cand) {fNumCandidate = cand;} - // virtual Int_t SetNumCandidateCut(Int_t candcut) {fNumCandidateCut = candcut;} - // Others virtual void OpenInputFiles() {} @@ -79,23 +76,32 @@ class AliJetReader : public TObject protected: AliJetReader(const AliJetReader& rJetReader); AliJetReader& operator = (const AliJetReader& rhsr); - TChain *fChain; // chain for reconstructed tracks - TClonesArray *fMomentumArray; // array of particle momenta - TClonesArray *fArrayMC; //! array of mc particles - TTask *fFillUnitArray; //! task list for filling the UnitArray - AliESDEvent *fESD; // pointer to esd - AliJetReaderHeader *fReaderHeader; // pointer to header - TArrayI fSignalFlag; // to flag if a particle comes from pythia or - // from the underlying event - TArrayI fCutFlag; // to flag if a particle passed the pt cut or not - TClonesArray *fUnitArray; // array of digit position and energy - TRefArray *fRefArray; // ! array of digit position and energy - TClonesArray *fUnitArrayNoCuts; // array of digit position and energy - Bool_t fArrayInitialised; // To check that array of units is initialised - AliJetFillUnitArrayTracks *fFillUAFromTracks; - AliJetFillUnitArrayEMCalDigits *fFillUAFromEMCalDigits; - Int_t fNumCandidate; - Int_t fNumCandidateCut; + + TChain *fChain; // chain for reconstructed tracks + TChain *fTree; // tree for reconstructed tracks + TClonesArray *fMomentumArray; // array of particle momenta + TClonesArray *fArrayMC; //! array of mc particles + TTask *fFillUnitArray; //! task list for filling the UnitArray + AliESDEvent *fESD; // pointer to esd + AliJetReaderHeader *fReaderHeader; // pointer to header + AliHeader *fAliHeader; // AliHeader + TArrayI fSignalFlag; // to flag if a particle comes from pythia or + // from the underlying event + TArrayI fCutFlag; // to flag if a particle passed the pt cut or not + TClonesArray *fUnitArray; // array of digit position and energy + TClonesArray *fUnitArrayNoCuts; //|| array of digit position and energy + Bool_t fArrayInitialised; // To check that array of units is initialised + AliJetFillUnitArrayTracks *fFillUAFromTracks; // For charged particle task + AliJetFillUnitArrayEMCalDigits *fFillUAFromEMCalDigits; // For neutral particle task + Int_t fNumCandidate; // Number of entries different from zero in unitarray + Int_t fNumCandidateCut; // Number of entries different from zero in unitarray + // which pass pt cut + AliJetHadronCorrectionv1 *fHadronCorrector; //! Pointer to hadronic correction + Int_t fHCorrection; // Hadron correction flag + Int_t fECorrection; // Electron correction flag + Bool_t fEFlag; // Electron correction flag + + ClassDef(AliJetReader,1) }; diff --git a/JETAN/AliJetReaderHeader.cxx b/JETAN/AliJetReaderHeader.cxx index 7d8d74c191a..3d687b71ffa 100755 --- a/JETAN/AliJetReaderHeader.cxx +++ b/JETAN/AliJetReaderHeader.cxx @@ -32,6 +32,7 @@ AliJetReaderHeader::AliJetReaderHeader(): fFirst(0), fLast(-1), fOption(0), + fCluster(0), fDebug(0), fDZ(0), fSignalPerBg(0), @@ -55,6 +56,7 @@ AliJetReaderHeader::AliJetReaderHeader(const char * name): fFirst(0), fLast(-1), fOption(0), + fCluster(0), fDebug(0), fDZ(0), fSignalPerBg(0), diff --git a/JETAN/AliJetReaderHeader.h b/JETAN/AliJetReaderHeader.h index cbf8208d1ed..38258d7e474 100755 --- a/JETAN/AliJetReaderHeader.h +++ b/JETAN/AliJetReaderHeader.h @@ -9,7 +9,7 @@ // // Author: jgcn@mda.cinvestav.mx //--------------------------------------------------------------------- - +#include #include #include @@ -35,8 +35,9 @@ class AliJetReaderHeader : public TNamed Int_t GetFirstEvent() const {return fFirst;} Int_t GetLastEvent() const {return fLast;} Int_t GetDetector() const {return fOption;} + Int_t GetCluster() const {return fCluster;} Bool_t GetDZ() const {return fDZ;} - Int_t GetDebug() const {return fDebug;} + Int_t GetDebug() const {cout << "coucou" << endl; return fDebug;} Int_t GetSignalPerBg() const {return fSignalPerBg;} // Setters @@ -54,12 +55,14 @@ class AliJetReaderHeader : public TNamed virtual void SetPtCut(Float_t par = 2.0) {fPtCut = par;} virtual void SetDZ(Bool_t deadzone = 0) {fDZ = deadzone;} virtual void SetDetector(Int_t option = 0) {fOption = option;} + virtual void SetCluster(Int_t option = 0) {fCluster = option;} virtual void SetDebug(Int_t debug = 0) {fDebug = debug;} protected: Int_t fFirst; // First and last events analyzed Int_t fLast; // in current set of files - Int_t fOption; // detector used for jet reconstruction + Int_t fOption; // detector used for jet reconstruction + Int_t fCluster; // cluster type Int_t fDebug; // debug option Bool_t fDZ; // include dead zones or not Int_t fSignalPerBg; diff --git a/JETAN/AliJetUnitArray.cxx b/JETAN/AliJetUnitArray.cxx index 2fea6528880..5b55049f268 100644 --- a/JETAN/AliJetUnitArray.cxx +++ b/JETAN/AliJetUnitArray.cxx @@ -13,14 +13,24 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -//------------------------------------------------------------------ + +/* $Id$ */ + +//_________________________________________________________________________ // Unit used by UA1 algorithm -// Authors: Sarah Blyth (LBL/UCT) -// Magali Estienne (IReS) (new version for JETAN) -//------------------------------------------------------------------ +// -- +//*-- Author: Sarah Blyth (LBL/UCT) +// -- +// Revised Version for JETAN +// -- Magali Estienne (IReS) -#include "AliJetUnitArray.h" +//#include + +#include +#include +#include +#include "AliJetUnitArray.h" ClassImp(AliJetUnitArray) @@ -36,16 +46,50 @@ AliJetUnitArray::AliJetUnitArray(): fUnitClusterID(0), fUnitFlag(kOutJet), fUnitCutFlag(kPtSmaller), + fUnitCutFlag2(kPtSmaller), fUnitSignalFlag(kBad), fUnitDetectorFlag(kTpc), fUnitPx(0.), fUnitPy(0.), fUnitPz(0.), - fUnitMass(0.) + fUnitMass(0.), + fV(0), + fVc(0), + fVn(0), + fVet(0) { // Default constructor } +AliJetUnitArray::AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t phi, Float_t en, Float_t Deta, Float_t Dphi, AliJetFinderUnitDetectorFlagType_t det, AliJetFinderUnitFlagType_t inout, AliJetFinderUnitCutFlagType_t cut, AliJetFinderUnitCutFlagType_t cut2, AliJetFinderUnitSignalFlagType_t signal,Float_t mass, Int_t clusId): + fUnitEnergy(en), + fUnitEta(eta), + fUnitPhi(phi), + fUnitDeta(Deta), + fUnitDphi(Dphi), + fUnitID(absId), + fUnitTrackID(esdId), + fUnitNum(0), + fUnitClusterID(clusId), + fUnitFlag(inout), + fUnitCutFlag(cut), + fUnitCutFlag2(cut2), + fUnitSignalFlag(signal), + fUnitDetectorFlag(det), + fUnitPx(0.), + fUnitPy(0.), + fUnitPz(0.), + fUnitMass(mass), + fV(0), + fVc(0), + fVn(0), + fVet(0) +{ + //abs ID (in a eta,phi grid, track ID in ESD, eta, phi, energy, px, py, pz, Deta, Dphi, detector flag, in/out jet, mass + + // Constructor 2 +} + AliJetUnitArray::AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t phi, Float_t en, Float_t px, Float_t py, Float_t pz, Float_t Deta, Float_t Dphi, AliJetFinderUnitDetectorFlagType_t det, AliJetFinderUnitFlagType_t inout, AliJetFinderUnitCutFlagType_t cut, Float_t mass, Int_t clusId): fUnitEnergy(en), fUnitEta(eta), @@ -58,22 +102,150 @@ AliJetUnitArray::AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t fUnitClusterID(clusId), fUnitFlag(inout), fUnitCutFlag(cut), - fUnitSignalFlag(kBad), + fUnitCutFlag2(kPtSmaller), + fUnitSignalFlag(kBad), fUnitDetectorFlag(det), fUnitPx(px), fUnitPy(py), fUnitPz(pz), - fUnitMass(mass) + fUnitMass(mass), + fV(0), + fVc(0), + fVn(0), + fVet(0) { // Constructor 2 } - + +//------------------------------------------------------------------------ AliJetUnitArray::~AliJetUnitArray() { // Destructor } - -Bool_t AliJetUnitArray::operator>(AliJetUnitArray &unit) const + +//------------------------------------------------------------------------ +void AliJetUnitArray::SetUnitSignalFlagC(Bool_t init, AliJetFinderUnitSignalFlagType_t flag) +{ + // Set signal flag of the charged particle + if(init){ + if(!fVc.empty()) + fVc.clear(); + } + else fVc.push_back(flag); +} + +//------------------------------------------------------------------------ +void AliJetUnitArray::SetUnitSignalFlagN(Bool_t init, AliJetFinderUnitSignalFlagType_t flag) +{ + // Set signal flag of the neutral cell + if(init){ + if(!fVn.empty()) + fVn.clear(); + } + else fVn.push_back(flag); +} + +//------------------------------------------------------------------------ +void AliJetUnitArray::SetUnitEtN(Bool_t init, Float_t et) +{ + // Set transverse energy of the neutral cell + if(init){ + if(!fVet.empty()) + fVet.clear(); + } + else fVet.push_back(et); +} + + +//------------------------------------------------------------------------ +void AliJetUnitArray::SetUnitPxPyPz(Bool_t init, vector v3) +{ + // Set momentum components of the charged particle + if(init) + { + if(!fV.empty()){ + fV.clear(); + } + } + else{ + fV.push_back(v3); + } +} + +//------------------------------------------------------------------------ +Bool_t AliJetUnitArray::GetUnitSignalFlagC(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagc) +{ + // Get signal flag of the charged particle + if(ind <= (Int_t)fVc.size()) + { + flagc = (AliJetFinderUnitSignalFlagType_t)fVc[ind]; + return kTRUE; + } + else return kFALSE; +} + +//------------------------------------------------------------------------ +Bool_t AliJetUnitArray::GetUnitSignalFlagN(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagn) +{ + // Get signal flag of the neutral cell + if(ind <= (Int_t)fVn.size()) + { + flagn = (AliJetFinderUnitSignalFlagType_t)fVn[ind]; + return kTRUE; + } + else return kFALSE; +} + +//------------------------------------------------------------------------ +Bool_t AliJetUnitArray::GetUnitEtN(Int_t ind, Float_t &et) +{ + // Get transverse energy of the neutral cell + if(ind <= (Int_t)fVet.size()) + { + et = (Float_t)fVet[ind]; + return kTRUE; + } + else return kFALSE; +} + +//------------------------------------------------------------------------ +Bool_t AliJetUnitArray::GetUnitPxPyPz(Int_t ind, Float_t &px, Float_t &py, Float_t &pz) +{ + // Get momentum components of the charged particle + if(ind <= (Int_t)fV.size()) + { + px = (Float_t)fV[ind][0]; + py = (Float_t)fV[ind][1]; + pz = (Float_t)fV[ind][2]; + return kTRUE; + } + else return kFALSE; +} + +//------------------------------------------------------------------------ +Bool_t AliJetUnitArray::GetUnitPxPyPzE(Int_t ind, Float_t &px, Float_t &py, Float_t &pz, Float_t &en) +{ +// Get 4-momentum components of the charged particle + if(ind <= (Int_t)fV.size()) + { + px = (Float_t)fV[ind][0]; + py = (Float_t)fV[ind][1]; + pz = (Float_t)fV[ind][2]; + en = TMath::Sqrt(px*px+py*py+pz*pz); + return kTRUE; + } + else return kFALSE; +} + +//------------------------------------------------------------------------ +Float_t AliJetUnitArray::EtaToTheta(Float_t arg) const +{ + // Eta to theta transformation + return 2.*atan(exp(-arg)); +} + +//------------------------------------------------------------------------ +Bool_t AliJetUnitArray::operator>(AliJetUnitArray unit) const { // Greater than operator used by sort if( fUnitEnergy > unit.GetUnitEnergy()) @@ -82,7 +254,8 @@ Bool_t AliJetUnitArray::operator>(AliJetUnitArray &unit) const return kFALSE; } -Bool_t AliJetUnitArray::operator<( AliJetUnitArray &unit) const +//------------------------------------------------------------------------ +Bool_t AliJetUnitArray::operator<( AliJetUnitArray unit) const { // Less than operator used by sort if( fUnitEnergy < unit.GetUnitEnergy()) @@ -91,7 +264,8 @@ Bool_t AliJetUnitArray::operator<( AliJetUnitArray &unit) const return kFALSE; } -Bool_t AliJetUnitArray::operator==( AliJetUnitArray &unit) const +//------------------------------------------------------------------------ +Bool_t AliJetUnitArray::operator==( AliJetUnitArray unit) const { // equality operator used by sort if( fUnitEnergy == unit.GetUnitEnergy()) diff --git a/JETAN/AliJetUnitArray.h b/JETAN/AliJetUnitArray.h index 858a1e70d7b..1673cda21be 100644 --- a/JETAN/AliJetUnitArray.h +++ b/JETAN/AliJetUnitArray.h @@ -4,27 +4,44 @@ /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * See cxx source for full Copyright notice */ -//------------------------------------------------------------------ -// Unit used by UA1 algorithm -// Authors: Sarah Blyth (LBL/UCT) -// Magali Estienne (IReS) (new version for JETAN) -//------------------------------------------------------------------ +/* $Id$ */ +// Class description : Unit red as input by jet finder algorithm to store +// the physical characteristics of a particle +// +// Author: magali.estienne@subatech.in2p3.fr +// +// Unit used by jet finder algorithm +// +// + +#include +#include #include +#include +#include #include "AliJetFinderTypes.h" class AliJetUnitArray : public TObject { public: AliJetUnitArray(); - AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t phi, Float_t en, Float_t px, Float_t py, Float_t pz, Float_t Deta, Float_t Dphi, AliJetFinderUnitDetectorFlagType_t det, AliJetFinderUnitFlagType_t inout, AliJetFinderUnitCutFlagType_t cut, Float_t mass, Int_t clusId); + AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t phi, Float_t en, Float_t Deta, Float_t Dphi, + AliJetFinderUnitDetectorFlagType_t det, AliJetFinderUnitFlagType_t inout, AliJetFinderUnitCutFlagType_t cut, + AliJetFinderUnitCutFlagType_t cut2, AliJetFinderUnitSignalFlagType_t signal, Float_t mass, Int_t clusId); + AliJetUnitArray(Int_t absId, Int_t esdId, Float_t eta, Float_t phi, Float_t en, Float_t px, Float_t py, Float_t pz, + Float_t Deta, Float_t Dphi, AliJetFinderUnitDetectorFlagType_t det, AliJetFinderUnitFlagType_t inout, + AliJetFinderUnitCutFlagType_t cut, Float_t mass, Int_t clusId); ~AliJetUnitArray(); // Setter void SetUnitEnergy(Float_t energy) {fUnitEnergy = energy;} void SetUnitEta(Float_t eta) {fUnitEta = eta;} void SetUnitPhi(Float_t phi) {fUnitPhi = phi;} + void SetUnitPx(Float_t px) {fUnitPx = px;} + void SetUnitPy(Float_t py) {fUnitPy = py;} + void SetUnitPz(Float_t pz) {fUnitPz = pz;} void SetUnitDeta(Float_t deta) {fUnitDeta = deta;} void SetUnitDphi(Float_t dphi) {fUnitDphi = dphi;} void SetUnitID(Int_t id) {fUnitID = id;} @@ -39,6 +56,10 @@ class AliJetUnitArray : public TObject { fUnitCutFlag = cutFlag; } + void SetUnitCutFlag2(AliJetFinderUnitCutFlagType_t cutFlag) + { + fUnitCutFlag2 = cutFlag; + } void SetUnitSignalFlag(AliJetFinderUnitSignalFlagType_t signalFlag) { fUnitSignalFlag = signalFlag; @@ -47,21 +68,36 @@ class AliJetUnitArray : public TObject { fUnitDetectorFlag = detectorflag; } - void SetUnitPxPyPz(Double_t *pxyz) {fUnitPx = pxyz[0]; fUnitPy = pxyz[1]; fUnitPz = pxyz[2];} + void SetUnitPxPyPz(Bool_t init, vector v3); + void SetUnitEtN(Bool_t init, Float_t et); // Added for background studies + void SetUnitSignalFlagC(Bool_t init, AliJetFinderUnitSignalFlagType_t flag); + void SetUnitSignalFlagN(Bool_t init, AliJetFinderUnitSignalFlagType_t flag); void SetUnitMass(Float_t mass) {fUnitMass = mass;} // Getter Float_t GetUnitEnergy() const {return fUnitEnergy;} Float_t GetUnitEta() const {return fUnitEta;} - Float_t GetUnitPhi() const {return fUnitPhi;} + Float_t GetUnitPhi() const {return fUnitPhi;} + Float_t GetUnitPx() const {return fUnitPx;} + Float_t GetUnitPy() const {return fUnitPy;} + Float_t GetUnitPz() const {return fUnitPz;} Float_t GetUnitDeta() const {return fUnitDeta;} Float_t GetUnitDphi() const {return fUnitDphi;} Int_t GetUnitID() const {return fUnitID;} - Int_t GetUnitTrackID() const {return fUnitTrackID;} + Int_t GetUnitTrackID() const {return fUnitTrackID;} Int_t GetUnitEntries() const {return fUnitNum;} Int_t GetUnitClusterID() const {return fUnitClusterID;} Float_t GetUnitMass() const {return fUnitMass;} - Bool_t GetUnitPxPyPz(Double_t* p) const {p[0]=fUnitPx; p[1]=fUnitPy; p[2]=fUnitPz; return kTRUE;} + Bool_t GetUnitPxPyPz(Int_t ind, Float_t &px, Float_t &py, Float_t &pz); + Bool_t GetUnitPxPyPzE(Int_t ind, Float_t &px, Float_t &py, Float_t &pz, Float_t &en); + Bool_t GetUnitEtN(Int_t ind, Float_t &et); // Added for background studies + Bool_t GetUnitSignalFlagC(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagc); + Bool_t GetUnitSignalFlagN(Int_t ind, AliJetFinderUnitSignalFlagType_t &flagn); + Int_t GetUnitVectorSize() {return fV.size();} + Int_t GetUnitVectorSizeC() {return fVc.size();} + Int_t GetUnitVectorSizeN() {return fVn.size();} + + Float_t EtaToTheta(Float_t arg) const; AliJetFinderUnitFlagType_t GetUnitFlag() const { @@ -71,6 +107,10 @@ class AliJetUnitArray : public TObject { return fUnitCutFlag; } + AliJetFinderUnitCutFlagType_t GetUnitCutFlag2() const + { + return fUnitCutFlag2; + } AliJetFinderUnitSignalFlagType_t GetUnitSignalFlag() const { return fUnitSignalFlag; @@ -80,35 +120,36 @@ class AliJetUnitArray : public TObject return fUnitDetectorFlag; } - - Bool_t operator> ( AliJetUnitArray &unit1) const; - Bool_t operator< ( AliJetUnitArray &unit1) const; - Bool_t operator== ( AliJetUnitArray &unit1) const; - protected: - Float_t fUnitEnergy; // Energy of the unit - Float_t fUnitEta; // Eta of the unit - Float_t fUnitPhi; // Phi of the unit - Float_t fUnitDeta; // Delta Eta of the unit - Float_t fUnitDphi; // Delta Phi of the unit - Int_t fUnitID; // ID of the unit - Int_t fUnitTrackID; // ID of the unit - Int_t fUnitNum; // number of units - Int_t fUnitClusterID; // ID of the unit + Bool_t operator> ( AliJetUnitArray unit1) const; + Bool_t operator< ( AliJetUnitArray unit1) const; + Bool_t operator== ( AliJetUnitArray unit1) const; + + Float_t fUnitEnergy; // Energy (Pt,et) of the unit + Float_t fUnitEta; // Eta of the unit + Float_t fUnitPhi; // Phi of the unit + Float_t fUnitDeta; // Delta Eta of the unit + Float_t fUnitDphi; // Delta Phi of the unit + Int_t fUnitID; // ID of the unit + Int_t fUnitTrackID; // ID of a given charged track + Int_t fUnitNum; // Number of units + Int_t fUnitClusterID; // ID for clusters AliJetFinderUnitFlagType_t fUnitFlag; // Flag of the unit - AliJetFinderUnitCutFlagType_t fUnitCutFlag; // Flag of the unit - AliJetFinderUnitSignalFlagType_t fUnitSignalFlag; // Flag of the unit + AliJetFinderUnitCutFlagType_t fUnitCutFlag; // Cut flag of the unit in the tpc + AliJetFinderUnitCutFlagType_t fUnitCutFlag2; // Cut flag of the unit in the emcal + AliJetFinderUnitSignalFlagType_t fUnitSignalFlag; // Signal flag of the unit AliJetFinderUnitDetectorFlagType_t fUnitDetectorFlag; // Detector flag of the unit - Float_t fUnitPx; // Px of charged track - Float_t fUnitPy; // Py of charged track - Float_t fUnitPz; // Pz of charged track - Float_t fUnitMass; // Mass of charged particle - - private: - AliJetUnitArray(const AliJetUnitArray &det); - AliJetUnitArray &operator=(const AliJetUnitArray &det); + Float_t fUnitPx; // Px of charged track + Float_t fUnitPy; // Py of charged track + Float_t fUnitPz; // Pz of charged track + Float_t fUnitMass; // Mass of particle + vector< vector< Float_t > > fV; //|| vector to store part information in each cell + vector< AliJetFinderUnitSignalFlagType_t > fVc; //|| added for background studies + vector< AliJetFinderUnitSignalFlagType_t > fVn; //|| added for background studies + vector< Float_t > fVet; //|| added for background studies ClassDef(AliJetUnitArray,1) + }; #endif diff --git a/JETAN/AliSISConeJetFinder.cxx b/JETAN/AliSISConeJetFinder.cxx index 307a2b77300..591a449c339 100644 --- a/JETAN/AliSISConeJetFinder.cxx +++ b/JETAN/AliSISConeJetFinder.cxx @@ -15,26 +15,31 @@ //--------------------------------------------------------------------- -// SISCone (FastJet v2.3.4) finder algorithm interface +// FastJet v2.3.4 finder algorithm interface // // Author: swensy.jangal@ires.in2p3.fr -// +// +// Last modification: Neutral cell energy included in the jet reconstruction +// +// Author: Magali.estienne@subatech.in2p3.fr //--------------------------------------------------------------------- + #include #include -#include #include #include #include #include #include +#include #include "AliHeader.h" #include "AliJet.h" #include "AliJetKineReader.h" #include "AliJetReader.h" #include "AliJetReaderHeader.h" +#include "AliJetUnitArray.h" #include "AliSISConeJetFinder.h" #include "AliSISConeJetHeader.h" @@ -43,21 +48,23 @@ #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" -// Get info on how fastjet was configured +// get info on how fastjet was configured #include "fastjet/config.h" #ifdef ENABLE_PLUGIN_SISCONE #include "fastjet/SISConePlugin.hh" #endif -#include // Needed for internal io +#include // needed for internal io #include #include using namespace std; + ClassImp(AliSISConeJetFinder) + //____________________________________________________________________________ AliSISConeJetFinder::AliSISConeJetFinder(): @@ -77,10 +84,10 @@ AliSISConeJetFinder::~AliSISConeJetFinder() void AliSISConeJetFinder::FindJets() { - Bool_t debug = kFALSE; - // Pick up siscone header AliSISConeJetHeader *header = (AliSISConeJetHeader*)fHeader; + Bool_t debug = header->GetDebug(); // debug option + Int_t fOpt = fReader->GetReaderHeader()->GetDetector(); // Check if we are reading AOD jets TRefArray *refs = 0; @@ -97,83 +104,252 @@ void AliSISConeJetFinder::FindJets() Double_t ptProtoJetMin = header->GetPtProtojetMin(); // pT min of protojets Double_t caching = header->GetCaching(); // do we record found cones for this set of data? - if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale +// if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale +// Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets fastjet::JetDefinition::Plugin * plugin; plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching); + vector inputParticles; //******************************** READING OF INPUT PARTICLES // Here we look for px, py pz and energy of each particle that we gather in a PseudoJet object, and we put all these PseudoJet in a vector of PseudoJets : input_particles. - TClonesArray *lvArray = fReader->GetMomentumArray(); - Int_t nIn = lvArray->GetEntries(); - - // We check if lvArray is ok - if(lvArray == 0) + if(fOpt==0) + { + TClonesArray *lvArray = fReader->GetMomentumArray(); + Int_t nIn = lvArray->GetEntries(); + + // We check if lvArray is ok + if(lvArray == 0) + { + cout << "Could not get the momentum array" << endl; + return; + } + + if(nIn == 0)// nIn = Number of particles in the event + { + if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; + return; + } + + fJets->SetNinput(nIn) ; // fJets = AliJet number of input objects + Float_t px,py,pz,en; + + // Load input vectors + for(Int_t i = 0; i < nIn; i++) + { + TLorentzVector *lv = (TLorentzVector*) lvArray->At(i); + px = lv->Px(); + py = lv->Py(); + pz = lv->Pz(); + en = lv->Energy(); + + fastjet::PseudoJet inputPart(px,py,pz,en); + inputPart.set_user_index(i); + inputParticles.push_back(inputPart); + + } + } + else { + TClonesArray* fUnit = fReader->GetUnitArray(); + if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; } + Int_t nCandidate = fReader->GetNumCandidate(); + Int_t nIn = fUnit->GetEntries(); + if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; } + fJets->SetNinput(nCandidate); // number of input objects // ME + // Information extracted from fUnitArray + // load input vectors and calculate total energy in array + Float_t pt,eta,phi,theta,px,py,pz,en; + Int_t ipart = 0; + for(Int_t i=0; iAt(i); + + if(uArray->GetUnitEnergy()>0.){ + + // It is not necessary anymore to cut on particle pt + pt = uArray->GetUnitEnergy(); + eta = uArray->GetUnitEta(); + phi = uArray->GetUnitPhi(); + theta = EtaToTheta(eta); + en = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta)); + px = TMath::Cos(phi)*pt; + py = TMath::Sin(phi)*pt; + pz = en*TMath::TanH(eta); + if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl; + + fastjet::PseudoJet inputPart(px,py,pz,en); // create PseudoJet object + inputPart.set_user_index(ipart); //label the particle into Fastjet algortihm + inputParticles.push_back(inputPart); // back of the input_particles vector + ipart++; + } + } // End loop on UnitArray + } + +//******************************** CHOICE OF JET AREA +// Here we determine jets area for subtracting background later +// For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez + + Double_t ghostEtamax = header->GetGhostEtaMax(); // maximum eta in which a ghost can be generated + Double_t ghostArea = header->GetGhostArea(); // area of a ghost + Int_t activeAreaRepeats = header->GetActiveAreaRepeats(); // do we repeat area calculation? + Double_t gridScatter = header->GetGridScatter(); // fractional random fluctuations of the position of the ghosts on the y-phi grid + Double_t ktScatter = header->GetKtScatter(); // fractional random fluctuations of the tranverse momentum of the ghosts on the y-phi grid + Double_t meanGhostKt = header->GetMeanGhostKt(); // average transverse momentum of the ghosts. + + Double_t areaTypeNumber = header->GetAreaTypeNumber(); // the number determines jet area type + fastjet::AreaType areaType = fastjet::active_area; + if (areaTypeNumber == 1) areaType = fastjet::active_area; + if (areaTypeNumber == 2) areaType = fastjet::active_area_explicit_ghosts; + if (areaTypeNumber == 3) areaType = fastjet::one_ghost_passive_area; + if (areaTypeNumber == 4) areaType = fastjet::passive_area; + if (areaTypeNumber == 5) areaType = fastjet::voronoi_area; + + fastjet::AreaDefinition areaDef; + + if (areaTypeNumber < 5) { - cout << "Could not get the momentum array" << endl; - return; + fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt); + areaDef = fastjet::AreaDefinition(areaType,ghostSpec); } - if(nIn == 0)// nIn = Number of particles in the event + if (areaTypeNumber == 5) { - if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; - return; + Double_t effectiveRFact = header->GetEffectiveRFact(); + fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact); + areaDef = fastjet::AreaDefinition(areaType,ghostSpec); } - Int_t nJets = 0; // Number of jets in this event - fJets->SetNinput(nIn) ; // fJets = AliJet number of input objects - Float_t px,py,pz,en; - vector inputParticles; +//******************************** +//******************************** - // Load input vectors - for(Int_t i = 0; i < nIn; i++) - { - TLorentzVector *lv = (TLorentzVector*) lvArray->At(i); - px = lv->Px(); - py = lv->Py(); - pz = lv->Pz(); - en = lv->Energy(); - - fastjet::PseudoJet inputPart(px,py,pz,en); - inputPart.set_user_index(i); - inputParticles.push_back(inputPart); - - } - -//******************************** JETS FINDING + Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not - fastjet::ClusterSequence clustSeq(inputParticles, plugin); +//******************************** +//******************************** -//***************************** JETS EXTRACTION AND CORRECTION +//******************************** +//******************************** BG SUBTRACTION + if (bgMode == 1)// BG subtraction + { + //******************************** JETS FINDING AND EXTRACTION + fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef); + // Here we extract inclusive jets with pt > ptmin, sorted by pt + Double_t ptMin = header->GetMinJetPt(); + vector inclusiveJets = clustSeq.inclusive_jets(ptMin); + vector jets = sorted_by_pt(inclusiveJets); + + //***************************** BACKGROUND SUBTRACTION + + // Set the rapidity-azimuth range within which to study background + Double_t rapMin = header->GetRapMin(); + Double_t rapMax = header->GetRapMax(); + Double_t phiMin = header->GetPhiMin(); + Double_t phiMax = header->GetPhiMax(); + fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax); + + // rho calculus from Cambridge/Aachen clustering (not from SISCONE as it gives too small areas) + Int_t algo = header->GetBGAlgorithm(); + fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm; + if (algo == 0) algorithm = fastjet::kt_algorithm; + if (algo == 1) algorithm = fastjet::cambridge_algorithm; + fastjet::JetDefinition jetDefForRho(algorithm, 0.5); + fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefForRho, areaDef); + Double_t rho = csForRho.median_pt_per_unit_area_4vector(range); + cout<<"rho = "< corrJets = clustSeq.subtracted_jets(rho, ptMin); + + //***************************** JETS DISPLAY + + for (size_t j = 0; j < jets.size(); j++) + { + // If the jet is only made of ghosts, continue. + if (clustSeq.is_pure_ghost(jets[j]) == 1) continue; + + // If the correction is > jet energy px = py = pz = e = 0 + if (corrJets[j].px() == 0 && corrJets[j].py() == 0 && corrJets[j].pz() == 0 && corrJets[j].E() == 0) continue; + + cout<<"********************************** Reconstructed jet(s) (non corrected)"< ptmin, sorted by pt - Double_t ptMin = header->GetMinJetPt(); - vector inclusiveJets = clustSeq.inclusive_jets(ptMin); - vector jets = sorted_by_pt(inclusiveJets); +//******************************** +//******************************** NO BG SUBTRACTION - for (Int_t j = 0; j < jets.size(); j++) + if (bgMode == 0)// No BG subtraction { - cout<<"********************************** Reconstructed jet(s) (non corrected)"< ptmin, sorted by pt + Double_t ptMin = header->GetMinJetPt(); + vector inclusiveJets = clustSeq.inclusive_jets(ptMin); + vector jets = sorted_by_pt(inclusiveJets); + + //***************************** JETS DISPLAY + + for (size_t k = 0; k < jets.size(); k++) + { + cout<<"********************************** Reconstructed jet(s) (non corrected)"<CreateTasks(tree); + +} + + //____________________________________________________________________________ void AliSISConeJetFinder::WriteJHeaderToFile() { fHeader->Write(); } + +//____________________________________________________________________________ + diff --git a/JETAN/AliSISConeJetFinder.h b/JETAN/AliSISConeJetFinder.h index 974b6d75965..5604237692b 100644 --- a/JETAN/AliSISConeJetFinder.h +++ b/JETAN/AliSISConeJetFinder.h @@ -32,7 +32,6 @@ #include "AliFastJetHeader.h" #include "AliJetFinder.h" - using namespace std; class AliSISConeJetFinder : public AliJetFinder @@ -42,12 +41,14 @@ class AliSISConeJetFinder : public AliJetFinder AliSISConeJetFinder(); ~AliSISConeJetFinder(); - void FindJets(); + void FindJets(); // others - void WriteJHeaderToFile(); - + void WriteJHeaderToFile(); + Float_t EtaToTheta(Float_t arg); + void InitTask(TChain* tree); + protected: AliSISConeJetFinder(const AliSISConeJetFinder& rfj); AliSISConeJetFinder& operator = (const AliSISConeJetFinder& rsfj); diff --git a/JETAN/AliSISConeJetHeader.cxx b/JETAN/AliSISConeJetHeader.cxx index 6e3ab6601d9..ecf1e2a29f4 100644 --- a/JETAN/AliSISConeJetHeader.cxx +++ b/JETAN/AliSISConeJetHeader.cxx @@ -23,8 +23,8 @@ #include #include -#include "fastjet/AreaDefinition.hh" #include "fastjet/ClusterSequenceArea.hh" +#include "fastjet/AreaDefinition.hh" #include "fastjet/JetDefinition.hh" #include "AliSISConeJetHeader.h" @@ -36,24 +36,30 @@ ClassImp(AliSISConeJetHeader) AliSISConeJetHeader::AliSISConeJetHeader(): AliJetHeader("AliSISConeJetHeader"), fActiveAreaRepeats(1), + fAreaTypeNumber(4), + fBGAlgo(1), + fBGMode(1), fCaching(0), - fConeRadius(0.4), + fConeRadius(0.7), + fDebug(0), fEffectiveRFact(1), + fGhostEtaMax(4.0), fGhostArea(0.05), - fGhostEtaMax(2.0), - fMinJetPt(0), + fGridScatter(1), + fKtScatter(0.1), + fMeanGhostKt(1e-100), + fMinJetPt(2), fNPassMax(0), - fOverlapThreshold(0.5), + fOverlapThreshold(0.75), fPhiMax(TMath::TwoPi()), fPhiMin(0), - fRapMax(-0.9), - fRapMin(0.9), fPtProtoJetMin(2), + fRapMax(0.9), + fRapMin(-0.9), fSplitMergeScaleNumber(0), - fSplitMergeStoppingScale(0) - + fSplitMergeStoppingScale(0) { - // Constructor + // Constructor } //////////////////////////////////////////////////////////////////////// @@ -71,14 +77,23 @@ void AliSISConeJetHeader::PrintParameters() const cout<<"Do we record cones of these events ? (0 = no, 1 = yes) = "< +#include #include #include @@ -45,65 +46,314 @@ ClassImp(AliUA1JetFinderV2) -//////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////// AliUA1JetFinderV2::AliUA1JetFinderV2() : - AliJetFinder(), - fLego(0), - fDebug(0), - fOpt(0) + AliJetFinder(), + fLego(0), + fDebug(0), + fOpt(0) { + // // Constructor + // } //////////////////////////////////////////////////////////////////////// - AliUA1JetFinderV2::~AliUA1JetFinderV2() - { - // destructor + // + // Destructor + // } //////////////////////////////////////////////////////////////////////// +void AliUA1JetFinderV2::FindJetsC() +{ + // + // Used to find jets using charged particle momentum information + // + // 1) Fill cell map array + // 2) calculate total energy and fluctuation level + // 3) Run algorithm + // 3.1) look centroides in cell map + // 3.2) calculate total energy in cones + // 3.3) flag as a possible jet + // 3.4) reorder cones by energy + // 4) subtract backg in accepted jets + // 5) fill AliJet list + + // Transform input to pt,eta,phi plus lego + + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + TClonesArray* lvArray = fReader->GetMomentumArray(); + Int_t nIn = lvArray->GetEntries(); + + if (nIn == 0) return; + + // local arrays for input + Float_t* ptT = new Float_t[nIn]; + Float_t* etaT = new Float_t[nIn]; + Float_t* phiT = new Float_t[nIn]; + Float_t* cFlagT = new Float_t[nIn]; // Temporarily added + Float_t* sFlagT = new Float_t[nIn]; // Temporarily added + Int_t* injet = new Int_t[nIn]; + + //total energy in array + Float_t etbgTotal = 0.0; + TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0); + + // load input vectors and calculate total energy in array + for (Int_t i = 0; i < nIn; i++){ + TLorentzVector *lv = (TLorentzVector*) lvArray->At(i); + ptT[i] = lv->Pt(); + etaT[i] = lv->Eta(); + phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi()); + cFlagT[i] = fReader->GetCutFlag(i); // Temporarily added + sFlagT[i] = fReader->GetSignalFlag(i); // Temporarily added peut-etre a mettre apres cut en pt !!! + + if (fReader->GetCutFlag(i) != 1) continue; + fLego->Fill(etaT[i], phiT[i], ptT[i]); + hPtTotal->Fill(ptT[i]); + etbgTotal+= ptT[i]; + } + + fJets->SetNinput(nIn); + + // calculate total energy and fluctuation in map + Double_t meanpt = hPtTotal->GetMean(); + Double_t ptRMS = hPtTotal->GetRMS(); + Double_t npart = hPtTotal->GetEntries(); + Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS); + + // arrays to hold jets + Float_t* etaJet = new Float_t[30]; + Float_t* phiJet = new Float_t[30]; + Float_t* etJet = new Float_t[30]; + Float_t* etsigJet = new Float_t[30]; //signal et in jet + Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable) + Int_t* ncellsJet = new Int_t[30]; + Int_t* multJet = new Int_t[30]; + //--- Added for jet reordering at the end of the jet finding procedure + Float_t* etaJetOk = new Float_t[30]; + Float_t* phiJetOk = new Float_t[30]; + Float_t* etJetOk = new Float_t[30]; + Float_t* etsigJetOk = new Float_t[30]; //signal et in jet + Float_t* etallJetOk = new Float_t[30]; // total et in jet (tmp variable) + Int_t* ncellsJetOk = new Int_t[30]; + Int_t* multJetOk = new Int_t[30]; + //-------------------------- + Int_t nJets; // to hold number of jets found by algorithm + Int_t nj; // number of jets accepted + Float_t prec = header->GetPrecBg(); + Float_t bgprec = 1; + while(bgprec > prec){ + //reset jet arrays in memory + memset(etaJet,0,sizeof(Float_t)*30); + memset(phiJet,0,sizeof(Float_t)*30); + memset(etJet,0,sizeof(Float_t)*30); + memset(etallJet,0,sizeof(Float_t)*30); + memset(etsigJet,0,sizeof(Float_t)*30); + memset(ncellsJet,0,sizeof(Int_t)*30); + memset(multJet,0,sizeof(Int_t)*30); + //--- Added for jet reordering at the end of the jet finding procedure + memset(etaJetOk,0,sizeof(Float_t)*30); + memset(phiJetOk,0,sizeof(Float_t)*30); + memset(etJetOk,0,sizeof(Float_t)*30); + memset(etallJetOk,0,sizeof(Float_t)*30); + memset(etsigJetOk,0,sizeof(Float_t)*30); + memset(ncellsJetOk,0,sizeof(Int_t)*30); + memset(multJetOk,0,sizeof(Int_t)*30); + //-------------------------- + nJets = 0; + nj = 0; + + // reset particles-jet array in memory + memset(injet,-1,sizeof(Int_t)*nIn); + //run cone algorithm finder + RunAlgoritmC(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet); + + //run background subtraction + if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event + nj = header->GetNAcceptJets(); + else + nj = nJets; + //subtract background + Float_t etbgTotalN = 0.0; //new background + if(header->GetBackgMode() == 1) // standard + // SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(header->GetBackgMode() == 2) //cone + SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(header->GetBackgMode() == 3) //ratio + SubtractBackgRatio(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(header->GetBackgMode() == 4) //statistic + SubtractBackgStat(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + //calc precision + if(etbgTotalN != 0.0) + bgprec = (etbgTotal - etbgTotalN)/etbgTotalN; + else + bgprec = 0; + etbgTotal = etbgTotalN; // update with new background estimation + } //end while + + // add jets to list + Int_t* idxjets = new Int_t[nj]; + Int_t nselectj = 0; + printf("Found %d jets \n", nj); + // Reorder jets by et in cone + Int_t * idx = new Int_t[nJets]; + TMath::Sort(nJets, etJet, idx); + for(Int_t p = 0; p < nJets; p++){ + etaJetOk[p] = etaJet[idx[p]]; + phiJetOk[p] = phiJet[idx[p]]; + etJetOk[p] = etJet[idx[p]]; + etallJetOk[p] = etJet[idx[p]]; + ncellsJetOk[p] = ncellsJet[idx[p]]; + multJetOk[p] = multJet[idx[p]]; + } + + for(Int_t kj=0; kj (header->GetJetEtaMax())) || + (etaJetOk[kj] < (header->GetJetEtaMin())) || + (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin + Float_t px, py,pz,en; // convert to 4-vector + px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]); + py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]); + pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj]))); + en = TMath::Sqrt(px * px + py * py + pz * pz); + fJets->AddJet(px, py, pz, en); + + AliAODJet jet(px, py, pz, en); + jet.Print(""); + + AddJet(jet); + + idxjets[nselectj] = kj; + nselectj++; + } -void AliUA1JetFinderV2::FindJets() + //add signal percentage and total signal in AliJets for analysis tool + Float_t* percentage = new Float_t[nselectj]; + Int_t* ncells = new Int_t[nselectj]; + Int_t* mult = new Int_t[nselectj]; + for(Int_t i = 0; i< nselectj; i++) + { + percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]]; + ncells[i] = ncellsJetOk[idxjets[i]]; + mult[i] = multJetOk[idxjets[i]]; + } + + //add particle-injet relationship /// + for(Int_t bj = 0; bj < nIn; bj++) + { + if(injet[bj] == -1) continue; //background particle + Int_t bflag = 0; + for(Int_t ci = 0; ci< nselectj; ci++){ + if(injet[bj] == idxjets[ci]){ + injet[bj]= ci; + bflag++; + break; + } + } + if(bflag == 0) injet[bj] = -1; // set as background particle + } + + fJets->SetNCells(ncells); + fJets->SetPtFromSignal(percentage); + fJets->SetMultiplicities(mult); + fJets->SetInJet(injet); + fJets->SetEtaIn(etaT); + fJets->SetPhiIn(phiT); + fJets->SetPtIn(ptT); + fJets->SetEtAvg(etbgTotal/(4*(header->GetLegoEtaMax())*TMath::Pi())); + + //delete + delete[] ptT; + delete[] etaT; + delete[] phiT; + delete[] cFlagT; + delete[] sFlagT; + delete[] injet; + delete[] hPtTotal; + delete[] etaJet; + delete[] phiJet; + delete[] etJet; + delete[] etsigJet; + delete[] etallJet; + delete[] ncellsJet; + delete[] multJet; + delete[] idxjets; + delete[] percentage; + delete[] ncells; + delete[] mult; + //--- Added for jet reordering + delete etaJetOk; + delete phiJetOk; + delete etJetOk; + delete etsigJetOk; + delete etallJetOk; + delete ncellsJetOk; + delete multJetOk; + //-------------------------- + +} +//////////////////////////////////////////////////////////////////////// +void AliUA1JetFinderV2::FindJets() { - //1) Fill cell map array - //2) calculate total energy and fluctuation level - //3) Run algorithm - // 3.1) look centroides in cell map - // 3.2) calculate total energy in cones - // 3.3) flag as a possible jet - // 3.4) reorder cones by energy - //4) subtract backg in accepted jets - //5) fill AliJet list + // + // Used to find jets using charged particle momentum information + // & neutral energy from calo cells + // + // 1) Fill cell map array + // 2) calculate total energy and fluctuation level + // 3) Run algorithm + // 3.1) look centroides in cell map + // 3.2) calculate total energy in cones + // 3.3) flag as a possible jet + // 3.4) reorder cones by energy + // 4) subtract backg in accepted jets + // 5) fill AliJet list // transform input to pt,eta,phi plus lego - AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; - TClonesArray* fUnit = fReader->GetUnitArray(); - Int_t nCandidate = fReader->GetNumCandidate(); - Int_t nIn = fUnit->GetEntries(); + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + TClonesArray* fUnit = fReader->GetUnitArray(); + Int_t nCand = fReader->GetNumCandidate(); + Int_t nCandCut = fReader->GetNumCandidateCut(); + Int_t nIn = fUnit->GetEntries(); + Float_t fPtMin = fReader->GetReaderHeader()->GetPtCut(); if (nIn == 0) return; + Int_t nCandidateCut = 0; + Int_t nCandidate = 0; + + nCandidate = nCand; + nCandidateCut = nCandCut; + // local arrays for input No Cuts // Both pt < ptMin and pt > ptMin - Float_t* enT = new Float_t[nCandidate]; Float_t* ptT = new Float_t[nCandidate]; + Float_t* en2T = new Float_t[nCandidate]; + Float_t* pt2T = new Float_t[nCandidate]; + Int_t* detT = new Int_t[nCandidate]; Float_t* etaT = new Float_t[nCandidate]; Float_t* phiT = new Float_t[nCandidate]; - Float_t* detaT = new Float_t[nCandidate]; - Float_t* dphiT = new Float_t[nCandidate]; Float_t* cFlagT = new Float_t[nCandidate]; + Float_t* cFlag2T = new Float_t[nCandidate]; + Float_t* sFlagT = new Float_t[nCandidate]; Float_t* cClusterT = new Float_t[nCandidate]; - Float_t* idT = new Float_t[nCandidate]; - Int_t loop1 = 0; - Int_t* injet = new Int_t[nCandidate]; - Int_t* sflag = new Int_t[nCandidate]; - + Int_t* vectT = new Int_t[nCandidate]; + Int_t loop1 = 0; + Int_t* injet = new Int_t[nCandidate]; + Int_t* sflag = new Int_t[nCandidate]; + vector< vector > pxT; + vector< vector > pyT; + vector< vector > pzT; //total energy in array Float_t etbgTotal = 0.0; @@ -114,46 +364,167 @@ void AliUA1JetFinderV2::FindJets() Float_t *etaCell = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray Float_t *phiCell = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray Int_t *flagCell = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray - + Float_t *etCell2 = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray + Float_t *etaCell2 = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray + Float_t *phiCell2 = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray + Int_t *flagCell2 = new Int_t[nIn]; //! Cell phi - Extracted from UnitArray + // Information extracted from fUnitArray - // load input vectors and calculate total energy in array + // Load input vectors and calculate total energy in array for(Int_t i=0; iAt(i); - if(uArray->GetUnitCutFlag()==1){ - etCell[i] = uArray->GetUnitEnergy(); - if (etCell[i] > 0.0) etCell[i] -= header->GetMinCellEt(); - if (etCell[i] < 0.0) etCell[i] = 0.; - etaCell[i] = uArray->GetUnitEta(); - phiCell[i] = uArray->GetUnitPhi(); - flagCell[i] = 0; // default - } - else { - etCell[i] = 0.; - etaCell[i] = uArray->GetUnitEta(); - phiCell[i] = uArray->GetUnitPhi(); - flagCell[i] = 0; - } - + if(uArray->GetUnitEnergy()>0.){ + vector vtmpx; + vector vtmpy; + vector vtmpz; + for(Int_t j=0; jGetUnitVectorSize();j++) + { + Float_t x=0.; Float_t y=0.; Float_t z=0.; + uArray->GetUnitPxPyPz(j,x,y,z); + vtmpx.push_back(x); + vtmpy.push_back(y); + vtmpz.push_back(z); + } + pxT.push_back(vtmpx); + pyT.push_back(vtmpy); + pzT.push_back(vtmpz); + vtmpx.clear(); + vtmpy.clear(); + vtmpz.clear(); ptT[loop1] = uArray->GetUnitEnergy(); - enT[loop1] = uArray->GetUnitEnergy(); + detT[loop1] = uArray->GetUnitDetectorFlag(); etaT[loop1] = uArray->GetUnitEta(); phiT[loop1] = uArray->GetUnitPhi(); - detaT[loop1] = uArray->GetUnitDeta(); - dphiT[loop1] = uArray->GetUnitDphi(); - cFlagT[loop1]= uArray->GetUnitCutFlag(); - idT[loop1] = uArray->GetUnitID(); - if(cFlagT[loop1] == 1) { - hPtTotal->Fill(ptT[loop1]); - // fLego->Fill(etaT[i], phiT[i], ptT[i]); - etbgTotal+= ptT[loop1]; + cFlagT[loop1]= uArray->GetUnitCutFlag(); // pt cut tpc + cFlag2T[loop1]= uArray->GetUnitCutFlag2(); // pt cut emcal + sFlagT[loop1]= uArray->GetUnitSignalFlag(); + vectT[loop1] = uArray->GetUnitVectorSize(); + if(cFlagT[loop1] == 1 || cFlag2T[loop1] == 1) { + pt2T[loop1] = 0.; + en2T[loop1] = 0.; + if(detT[loop1]==1){ + en2T[loop1] = ptT[loop1] - header->GetMinCellEt(); + if(en2T[loop1] < 0) en2T[loop1]=0; + hPtTotal->Fill(en2T[loop1]); + etbgTotal += en2T[loop1]; + } + if(detT[loop1]==0){ // TPC+ITS + Float_t pt = 0.; + for(Int_t j=0; jGetUnitPxPyPz(j,x,y,z); + pt = TMath::Sqrt(x*x+y*y); + if(pt>fPtMin) { + pt2T[loop1] += pt; + en2T[loop1] += pt; + hPtTotal->Fill(pt); + etbgTotal+= pt; + } + } + } + if(detT[loop1]==2) { // EMCal + Float_t ptCTot = 0.; + Float_t pt = 0.; + Float_t enC = 0.; + for(Int_t j=0; jGetUnitPxPyPz(j,x,y,z); + pt = TMath::Sqrt(x*x+y*y); + if(pt>fPtMin) { + pt2T[loop1]+=pt; + en2T[loop1]+=pt; + hPtTotal->Fill(pt); + etbgTotal+= pt; + } + ptCTot += pt; + } + enC = ptT[loop1] - ptCTot - header->GetMinCellEt(); + if(enC < 0.) enC=0.; + en2T[loop1] += enC; + hPtTotal->Fill(enC); + etbgTotal+= enC; + } } loop1++; } - } - // fJets->SetNinput(nIn); + if(uArray->GetUnitCutFlag()==1) { + if(uArray->GetUnitDetectorFlag()==1){ // EMCal case + etCell[i] = uArray->GetUnitEnergy() - header->GetMinCellEt(); + if ((uArray->GetUnitEnergy() - header->GetMinCellEt()) < 0.0) etCell[i]=0.; + etaCell[i] = uArray->GetUnitEta(); + phiCell[i] = uArray->GetUnitPhi(); + flagCell[i] = 0; // default + etCell2[i] = etCell[i]; + etaCell2[i] = uArray->GetUnitEta(); + phiCell2[i] = uArray->GetUnitPhi(); + flagCell2[i] = 0; // default + } + if(uArray->GetUnitDetectorFlag()==0){ // TPC case + Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.; + for(Int_t j=0; jGetUnitVectorSize();j++) + { + Float_t x=0.; Float_t y=0.; Float_t z=0.; + uArray->GetUnitPxPyPz(j,x,y,z); + pt = TMath::Sqrt(x*x+y*y); + if(pt>fPtMin) { + et1 += pt; + et2 += pt; + } + } + etCell[i] = et1; + etCell2[i] = et2; + if(et1 < 0.) etCell[i] = etCell2[i] = 0.; + etaCell[i] = uArray->GetUnitEta(); + phiCell[i] = uArray->GetUnitPhi(); + flagCell[i] = 0; // default + etaCell2[i] = uArray->GetUnitEta(); + phiCell2[i] = uArray->GetUnitPhi(); + flagCell2[i] = 0; // default + } + if(uArray->GetUnitDetectorFlag()==2){ // TPC + EMCal case + Float_t ptCTot = 0.; + Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.; + Float_t enC = 0.; + for(Int_t j=0; jGetUnitVectorSize();j++) + { + Float_t x=0.; Float_t y=0.; Float_t z=0.; + uArray->GetUnitPxPyPz(j,x,y,z); + pt = TMath::Sqrt(x*x+y*y); + if(pt>fPtMin) { + et1 += pt; + et2 += pt; + } + ptCTot += pt; + } + enC = uArray->GetUnitEnergy() - ptCTot; + etCell[i] = et1 + enC - header->GetMinCellEt(); + etCell2[i] = et2 + enC - header->GetMinCellEt(); + if((enC + et1 - header->GetMinCellEt()) < 0.) etCell[i] = etCell2[i] = 0.; + etaCell[i] = uArray->GetUnitEta(); + phiCell[i] = uArray->GetUnitPhi(); + flagCell[i] = 0; // default + etaCell2[i] = uArray->GetUnitEta(); + phiCell2[i] = uArray->GetUnitPhi(); + flagCell2[i] = 0; // default + } + } + else { + etCell[i] = 0.; + etaCell[i] = uArray->GetUnitEta(); + phiCell[i] = uArray->GetUnitPhi(); + flagCell[i] = 0; + etCell2[i] = 0.; + etaCell2[i] = uArray->GetUnitEta(); + phiCell2[i] = uArray->GetUnitPhi(); + flagCell2[i] = 0; + } + } // end loop on nCandidate + fJets->SetNinput(nCandidate); // calculate total energy and fluctuation in map @@ -163,68 +534,110 @@ void AliUA1JetFinderV2::FindJets() Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS); // arrays to hold jets - Float_t* etaJet = new Float_t[30]; - Float_t* phiJet = new Float_t[30]; - Float_t* etJet = new Float_t[30]; + Float_t* etaJet = new Float_t[30]; + Float_t* phiJet = new Float_t[30]; + Float_t* etJet = new Float_t[30]; Float_t* etsigJet = new Float_t[30]; //signal et in jet - Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable) - Int_t* ncellsJet = new Int_t[30]; - Int_t* multJet = new Int_t[30]; - Int_t nJets; // to hold number of jets found by algorithm - Int_t nj; // number of jets accepted - Float_t prec = header->GetPrecBg(); - Float_t bgprec = 1; + Float_t* etallJet = new Float_t[30]; // total et in jet (tmp variable) + Int_t* ncellsJet = new Int_t[30]; + Int_t* multJet = new Int_t[30]; + //--- Added by me for jet reordering at the end of the jet finding procedure + Float_t* etaJetOk = new Float_t[30]; + Float_t* phiJetOk = new Float_t[30]; + Float_t* etJetOk = new Float_t[30]; + Float_t* etsigJetOk = new Float_t[30]; //signal et in jet + Float_t* etallJetOk = new Float_t[30]; // total et in jet (tmp variable) + Int_t* ncellsJetOk = new Int_t[30]; + Int_t* multJetOk = new Int_t[30]; + //-------------------------- + Int_t nJets; // to hold number of jets found by algorithm + Int_t nj; // number of jets accepted + Float_t prec = header->GetPrecBg(); + Float_t bgprec = 1; + while(bgprec > prec){ - //reset jet arrays in memory - memset(etaJet,0,sizeof(Float_t)*30); - memset(phiJet,0,sizeof(Float_t)*30); - memset(etJet,0,sizeof(Float_t)*30); - memset(etallJet,0,sizeof(Float_t)*30); - memset(etsigJet,0,sizeof(Float_t)*30); - memset(ncellsJet,0,sizeof(Int_t)*30); - memset(multJet,0,sizeof(Int_t)*30); - nJets = 0; - nj = 0; - // reset particles-jet array in memory - memset(injet,-1,sizeof(Int_t)*nCandidate); - //run cone algorithm finder - RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet); - //run background subtraction - if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event - nj = header->GetNAcceptJets(); - else - nj = nJets; - //subtract background - Float_t etbgTotalN = 0.0; //new background - if(header->GetBackgMode() == 1) // standar - SubtractBackg(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); - if(header->GetBackgMode() == 2) //cone - SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); - if(header->GetBackgMode() == 3) //ratio - SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); - if(header->GetBackgMode() == 4) //statistic - SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet); - //calc precision - if(etbgTotalN != 0.0) - bgprec = (etbgTotal - etbgTotalN)/etbgTotalN; - else - bgprec = 0; - etbgTotal = etbgTotalN; // update with new background estimation - } //end while + //reset jet arrays in memory + memset(etaJet,0,sizeof(Float_t)*30); + memset(phiJet,0,sizeof(Float_t)*30); + memset(etJet,0,sizeof(Float_t)*30); + memset(etallJet,0,sizeof(Float_t)*30); + memset(etsigJet,0,sizeof(Float_t)*30); + memset(ncellsJet,0,sizeof(Int_t)*30); + memset(multJet,0,sizeof(Int_t)*30); + //--- Added by me for jet reordering at the end of the jet finding procedure + memset(etaJetOk,0,sizeof(Float_t)*30); + memset(phiJetOk,0,sizeof(Float_t)*30); + memset(etJetOk,0,sizeof(Float_t)*30); + memset(etallJetOk,0,sizeof(Float_t)*30); + memset(etsigJetOk,0,sizeof(Float_t)*30); + memset(ncellsJetOk,0,sizeof(Int_t)*30); + memset(multJetOk,0,sizeof(Int_t)*30); + + nJets = 0; + nj = 0; + + // reset particles-jet array in memory + memset(injet,-1,sizeof(Int_t)*nCandidate); + //run cone algorithm finder + RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etCell2,etaCell2,phiCell2, + flagCell2,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet, + etallJet,ncellsJet); + + //run background subtraction + if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event + nj = header->GetNAcceptJets(); + else + nj = nJets; + + //subtract background + Float_t etbgTotalN = 0.0; //new background + if(header->GetBackgMode() == 1) // standard + SubtractBackg(nCandidate,nj,etbgTotalN,en2T,vectT,etaT,phiT,cFlagT,cFlag2T,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + // To be modified ------------------------ + if(header->GetBackgMode() == 2) //cone + SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(header->GetBackgMode() == 3) //ratio + SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + if(header->GetBackgMode() == 4) //statistic + SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet); + //---------------------------------------- + //calc precision + if(etbgTotalN != 0.0) + bgprec = (etbgTotal - etbgTotalN)/etbgTotalN; + else + bgprec = 0; + etbgTotal = etbgTotalN; // update with new background estimation + } //end while + // add jets to list Int_t* idxjets = new Int_t[nj]; Int_t nselectj = 0; printf("Found %d jets \n", nj); - for(Int_t kj=0; kj (header->GetJetEtaMax())) || - (etaJet[kj] < (header->GetJetEtaMin())) || - (etJet[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin + // Reorder jets by et in cone + // Sort jets by energy + Int_t * idx = new Int_t[nJets]; + TMath::Sort(nJets, etJet, idx); + for(Int_t p = 0; p < nJets; p++) + { + etaJetOk[p] = etaJet[idx[p]]; + phiJetOk[p] = phiJet[idx[p]]; + etJetOk[p] = etJet[idx[p]]; + etallJetOk[p] = etJet[idx[p]]; + ncellsJetOk[p] = ncellsJet[idx[p]]; + multJetOk[p] = multJet[idx[p]]; + } + + for(Int_t kj=0; kj (header->GetJetEtaMax())) || + (etaJetOk[kj] < (header->GetJetEtaMin())) || + (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin Float_t px, py,pz,en; // convert to 4-vector - px = etJet[kj] * TMath::Cos(phiJet[kj]); - py = etJet[kj] * TMath::Sin(phiJet[kj]); - pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj]))); + px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]); + py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]); + pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj]))); en = TMath::Sqrt(px * px + py * py + pz * pz); fJets->AddJet(px, py, pz, en); AliAODJet jet(px, py, pz, en); @@ -234,30 +647,34 @@ void AliUA1JetFinderV2::FindJets() idxjets[nselectj] = kj; nselectj++; - } + } + //add signal percentage and total signal in AliJets for analysis tool Float_t* percentage = new Float_t[nselectj]; Int_t* ncells = new Int_t[nselectj]; Int_t* mult = new Int_t[nselectj]; - for(Int_t i = 0; i< nselectj; i++){ - percentage[i] = etsigJet[idxjets[i]]/etJet[idxjets[i]]; - ncells[i] = ncellsJet[idxjets[i]]; - mult[i] = multJet[idxjets[i]]; + for(Int_t i = 0; i< nselectj; i++) + { + percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]]; + ncells[i] = ncellsJetOk[idxjets[i]]; + mult[i] = multJetOk[idxjets[i]]; + } + + //add particle-injet relationship /// + for(Int_t bj = 0; bj < nCandidate; bj++) + { + if(injet[bj] == -1) continue; //background particle + Int_t bflag = 0; + for(Int_t ci = 0; ci< nselectj; ci++){ + if(injet[bj] == idxjets[ci]){ + injet[bj]= ci; + bflag++; + break; + } + } + if(bflag == 0) injet[bj] = -1; // set as background particle } - //add particle-injet relationship /// - // for(Int_t bj = 0; bj < nIn; bj++){ - for(Int_t bj = 0; bj < nCandidate; bj++){ - if(injet[bj] == -1) continue; //background particle - Int_t bflag = 0; - for(Int_t ci = 0; ci< nselectj; ci++){ - if(injet[bj] == idxjets[ci]){ - injet[bj]= ci; - bflag++; - break; - } - } - if(bflag == 0) injet[bj] = -1; // set as background particle - } + fJets->SetNCells(ncells); fJets->SetPtFromSignal(percentage); fJets->SetMultiplicities(mult); @@ -265,19 +682,28 @@ void AliUA1JetFinderV2::FindJets() fJets->SetEtaIn(etaT); fJets->SetPhiIn(phiT); fJets->SetPtIn(ptT); - fJets->SetEtAvg(etbgTotal/(4*(header->GetLegoEtaMax())*TMath::Pi())); - + fJets->SetVectorSizeIn(vectT); + fJets->SetVectorPxIn(pxT); + fJets->SetVectorPyIn(pyT); + fJets->SetVectorPzIn(pzT); + fJets->SetDetectorFlagIn(detT); + fJets->SetEtAvg(etbgTotal/(2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax()-header->GetLegoPhiMin()))); //delete - delete enT; delete ptT; + delete en2T; + delete pt2T; delete etaT; delete phiT; - delete detaT; - delete dphiT; + pxT.clear(); + pyT.clear(); + pzT.clear(); + delete detT; delete cFlagT; + delete cFlag2T; + delete sFlagT; delete cClusterT; - delete idT; + delete vectT; delete injet; delete sflag; delete hPtTotal; @@ -285,6 +711,10 @@ void AliUA1JetFinderV2::FindJets() delete etaCell; delete phiCell; delete flagCell; + delete etCell2; + delete etaCell2; + delete phiCell2; + delete flagCell2; delete etaJet; delete phiJet; delete etJet; @@ -292,25 +722,41 @@ void AliUA1JetFinderV2::FindJets() delete etallJet; delete ncellsJet; delete multJet; + //--- Added for jet reordering + delete etaJetOk; + delete phiJetOk; + delete etJetOk; + delete etsigJetOk; + delete etallJetOk; + delete ncellsJetOk; + delete multJetOk; + //-------------------------- delete idxjets; delete percentage; delete ncells; delete mult; - } //////////////////////////////////////////////////////////////////////// - void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, - Int_t* flagCell, Float_t etbgTotal, Double_t dEtTotal, + Int_t* flagCell, Float_t* etCell2, Float_t* etaCell2, Float_t* phiCell2, + Int_t* flagCell2, Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etallJet, Int_t* ncellsJet) { + + Int_t nCell = nIn; - Int_t nCell = nIn; - + // Dump lego + // Check enough space! *to be done* AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + for(Int_t i=0; iGetMinMove(); @@ -318,19 +764,19 @@ void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell Float_t rc = header->GetRadius(); Float_t etseed = header->GetEtSeed(); - // tmp array of jets form algoritm + // Tmp array of jets form algoritm Float_t etaAlgoJet[30]; Float_t phiAlgoJet[30]; Float_t etAlgoJet[30]; Int_t ncellsAlgoJet[30]; - //run algorithm// + // Run algorithm// - // sort cells by et + // Sort cells by et Int_t * index = new Int_t[nCell]; TMath::Sort(nCell, etCell, index); - // variable used in centroide loop + // Variable used in centroide loop Float_t eta = 0.0; Float_t phi = 0.0; Float_t eta0 = 0.0; @@ -348,11 +794,12 @@ void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell Float_t phisb = 0.0; Float_t dphib = 0.0; - - for(Int_t icell = 0; icell < nCell; icell++){ + for(Int_t icell = 0; icell < nCell; icell++) + { Int_t jcell = index[icell]; if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed if(flagCell[jcell] != 0) continue; // if cell was used before + eta = etaCell[jcell]; phi = phiCell[jcell]; eta0 = eta; @@ -365,189 +812,604 @@ void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell etsb = ets; etasb = 0.0; phisb = 0.0; - for(Int_t kcell =0; kcell < nCell; kcell++){ + for(Int_t kcell =0; kcell < nCell; kcell++) + { Int_t lcell = index[kcell]; if(lcell == jcell) continue; // cell itself if(flagCell[lcell] != 0) continue; // cell used before - if(etCell[lcell] > etCell[jcell]) continue; + if(etCell[lcell] > etCell[jcell]) continue; // can this happen //calculate dr deta = etaCell[lcell] - eta; - dphi = phiCell[lcell] - phi; - if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); + dphi = TMath::Abs(phiCell[lcell] - phi); if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; dr = TMath::Sqrt(deta * deta + dphi * dphi); if(dr <= rc){ - // calculate offset from initiate cell - deta = etaCell[lcell] - eta0; - dphi = phiCell[lcell] - phi0; - if (dphi < -TMath::Pi()) dphi = dphi + 2.0 * TMath::Pi(); - if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi(); - etas = etas + etCell[lcell]*deta; - phis = phis + etCell[lcell]*dphi; - ets = ets + etCell[lcell]; - //new weighted eta and phi including this cell - eta = eta0 + etas/ets; - phi = phi0 + phis/ets; - // if cone does not move much, just go to next step - dphib = TMath::Abs(phi - phib); - if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib; - dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib); - if(dr <= minmove) break; - // cone should not move more than max_mov - dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets)); - if(dr > maxmove){ - eta = etab; - phi = phib; - ets = etsb; - etas = etasb; - phis = phisb; - }else{ // store this loop information - etab=eta; - phib=phi; - etsb = ets; - etasb = etas; - phisb = phis; - } - } + // calculate offset from initiate cell + deta = etaCell[lcell] - eta0; + dphi = phiCell[lcell] - phi0; + if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi(); + if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi(); + etas = etas + etCell[lcell]*deta; + phis = phis + etCell[lcell]*dphi; + ets = ets + etCell[lcell]; + //new weighted eta and phi including this cell + eta = eta0 + etas/ets; + phi = phi0 + phis/ets; + // if cone does not move much, just go to next step + dphib = TMath::Abs(phi - phib); + if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib; + dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib); + if(dr <= minmove) break; + // cone should not move more than max_mov + dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets)); + if(dr > maxmove){ + eta = etab; + phi = phib; + ets = etsb; + etas = etasb; + phis = phisb; + } else { // store this loop information + etab = eta; + phib = phi; + etsb = ets; + etasb = etas; + phisb = phis; + } + } // inside cone }//end of cells loop looking centroide //avoid cones overloap (to be implemented in the future) //flag cells in Rc, estimate total energy in cone - Float_t etCone = 0.0; - Int_t nCellIn = 0; - rc = header->GetRadius(); - for(Int_t ncell =0; ncell < nCell; ncell++){ - if(flagCell[ncell] != 0) continue; // cell used before - //calculate dr - deta = etaCell[ncell] - eta; - dphi = phiCell[ncell] - phi; - if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); - if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; - dr = TMath::Sqrt(deta * deta + dphi * dphi); - if(dr <= rc){ // cell in cone - flagCell[ncell] = -1; - etCone+=etCell[ncell]; - nCellIn++; - } + Float_t etCone = 0.0; + Int_t nCellIn = 0; + Int_t nCellOut = 0; + rc = header->GetRadius(); + + for(Int_t ncell =0; ncell < nCell; ncell++) + { + if(flagCell[ncell] != 0) continue; // cell used before + //calculate dr + deta = etaCell[ncell] - eta; + // if(deta <= rc){ // Added to improve velocity -> to be tested + dphi = phiCell[ncell] - phi; + if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); + if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + // if(dphi <= rc){ // Added to improve velocity -> to be tested + dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // cell in cone + flagCell[ncell] = -1; + etCone+=etCell[ncell]; + nCellIn++; + } + else nCellOut++; + // } // end deta <= rc + // } // end dphi <= rc } - // select jets with et > background - // estimate max fluctuation of background in cone - Double_t ncellin = (Double_t)nCellIn; - Double_t ntcell = (Double_t)nCell; - Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell); - // min cone et - Double_t etcmin = etCone ; // could be used etCone - etmin !! - //desicions !! etbmax < etcmin - for(Int_t mcell =0; mcell < nCell; mcell++){ - if(flagCell[mcell] == -1){ - if(etbmax < etcmin) - flagCell[mcell] = 1; //flag cell as used - else - flagCell[mcell] = 0; // leave it free - } - } - //store tmp jet info !!! - if(etbmax < etcmin) { - etaAlgoJet[nJets] = eta; - phiAlgoJet[nJets] = phi; - etAlgoJet[nJets] = etCone; - ncellsAlgoJet[nJets] = nCellIn; - nJets++; + // select jets with et > background + // estimate max fluctuation of background in cone + Double_t ncellin = (Double_t)nCellIn; + Double_t ntcell = (Double_t)nCell; + Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/(ntcell)); + // min cone et + Double_t etcmin = etCone ; // could be used etCone - etmin !! + //decisions !! etbmax < etcmin + + for(Int_t mcell =0; mcell < nCell; mcell++) + { + if(flagCell[mcell] == -1){ + if(etbmax < etcmin) + flagCell[mcell] = 1; //flag cell as used + else + flagCell[mcell] = 0; // leave it free + } } + //store tmp jet info !!! + if(etbmax < etcmin) + { + etaAlgoJet[nJets] = eta; + phiAlgoJet[nJets] = phi; + etAlgoJet[nJets] = etCone; + ncellsAlgoJet[nJets] = nCellIn; + nJets++; + } + + } // end of cells loop + + for(Int_t p = 0; p < nJets; p++) + { + etaJet[p] = etaAlgoJet[p]; + phiJet[p] = phiAlgoJet[p]; + etJet[p] = etAlgoJet[p]; + etallJet[p] = etAlgoJet[p]; + ncellsJet[p] = ncellsAlgoJet[p]; + } + + //delete + delete index; + +} + +//////////////////////////////////////////////////////////////////////// +void AliUA1JetFinderV2::RunAlgoritmC(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets, + Float_t* etJet,Float_t* etaJet, Float_t* phiJet, + Float_t* etallJet, Int_t* ncellsJet) +{ + // Dump lego + // Check enough space! *to be done* + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + Float_t etCell[60000]; //! Cell Energy + Float_t etaCell[60000]; //! Cell eta + Float_t phiCell[60000]; //! Cell phi + Int_t flagCell[60000]; //! Cell flag + + Int_t nCell = 0; + TAxis* xaxis = fLego->GetXaxis(); + TAxis* yaxis = fLego->GetYaxis(); + Float_t e = 0.0; + for (Int_t i = 1; i <= header->GetLegoNbinEta(); i++) + { + for (Int_t j = 1; j <= header->GetLegoNbinPhi(); j++) + { + e = fLego->GetBinContent(i,j); + if (e < 0.0) continue; // don't include this cells + Float_t eta = xaxis->GetBinCenter(i); + Float_t phi = yaxis->GetBinCenter(j); + etCell[nCell] = e; + etaCell[nCell] = eta; + phiCell[nCell] = phi; + flagCell[nCell] = 0; //default + nCell++; + } + } + + // Parameters from header + Float_t minmove = header->GetMinMove(); + Float_t maxmove = header->GetMaxMove(); + Float_t rc = header->GetRadius(); + Float_t etseed = header->GetEtSeed(); + + // Tmp array of jets form algoritm + Float_t etaAlgoJet[30]; + Float_t phiAlgoJet[30]; + Float_t etAlgoJet[30]; + Int_t ncellsAlgoJet[30]; + + // Run algorithm// + + // Sort cells by et + Int_t * index = new Int_t[nCell]; + TMath::Sort(nCell, etCell, index); + // variable used in centroide loop + Float_t eta = 0.0; + Float_t phi = 0.0; + Float_t eta0 = 0.0; + Float_t phi0 = 0.0; + Float_t etab = 0.0; + Float_t phib = 0.0; + Float_t etas = 0.0; + Float_t phis = 0.0; + Float_t ets = 0.0; + Float_t deta = 0.0; + Float_t dphi = 0.0; + Float_t dr = 0.0; + Float_t etsb = 0.0; + Float_t etasb = 0.0; + Float_t phisb = 0.0; + Float_t dphib = 0.0; - } // end of cells loop + for(Int_t icell = 0; icell < nCell; icell++) + { + Int_t jcell = index[icell]; + if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed + if(flagCell[jcell] != 0) continue; // if cell was used before + + eta = etaCell[jcell]; + phi = phiCell[jcell]; + eta0 = eta; + phi0 = phi; + etab = eta; + phib = phi; + ets = etCell[jcell]; + etas = 0.0; + phis = 0.0; + etsb = ets; + etasb = 0.0; + phisb = 0.0; + for(Int_t kcell =0; kcell < nCell; kcell++) + { + Int_t lcell = index[kcell]; + if(lcell == jcell) continue; // cell itself + if(flagCell[lcell] != 0) continue; // cell used before + if(etCell[lcell] > etCell[jcell]) continue; // can this happen + //calculate dr + deta = etaCell[lcell] - eta; + dphi = TMath::Abs(phiCell[lcell] - phi); + if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc) + { + // calculate offset from initiate cell + deta = etaCell[lcell] - eta0; + dphi = phiCell[lcell] - phi0; + if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi(); + if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi(); + etas = etas + etCell[lcell]*deta; + phis = phis + etCell[lcell]*dphi; + ets = ets + etCell[lcell]; + //new weighted eta and phi including this cell + eta = eta0 + etas/ets; + phi = phi0 + phis/ets; + // if cone does not move much, just go to next step + dphib = TMath::Abs(phi - phib); + if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib; + dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib); + if(dr <= minmove) break; + // cone should not move more than max_mov + dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets)); + if(dr > maxmove){ + eta = etab; + phi = phib; + ets = etsb; + etas = etasb; + phis = phisb; + } else { // store this loop information + etab=eta; + phib=phi; + etsb = ets; + etasb = etas; + phisb = phis; + } + } // inside cone + }//end of cells loop looking centroide + + // Avoid cones overloap (to be implemented in the future) + + // Flag cells in Rc, estimate total energy in cone + Float_t etCone = 0.0; + Int_t nCellIn = 0; + Int_t nCellOut = 0; + rc = header->GetRadius(); + for(Int_t ncell =0; ncell < nCell; ncell++) + { + if(flagCell[ncell] != 0) continue; // cell used before + //calculate dr + deta = etaCell[ncell] - eta; + dphi = phiCell[ncell] - phi; + if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi(); + if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // cell in cone + flagCell[ncell] = -1; + etCone+=etCell[ncell]; + nCellIn++; + } + else nCellOut++; + } + + // Select jets with et > background + // estimate max fluctuation of background in cone + Double_t ncellin = (Double_t)nCellIn; + Double_t ntcell = (Double_t)nCell; + Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell); + // min cone et + Double_t etcmin = etCone ; // could be used etCone - etmin !! + //decisions !! etbmax < etcmin + + for(Int_t mcell =0; mcell < nCell; mcell++){ + if(flagCell[mcell] == -1){ + if(etbmax < etcmin) + flagCell[mcell] = 1; //flag cell as used + else + flagCell[mcell] = 0; // leave it free + } + } + //store tmp jet info !!! + + if(etbmax < etcmin) { + etaAlgoJet[nJets] = eta; + phiAlgoJet[nJets] = phi; + etAlgoJet[nJets] = etCone; + ncellsAlgoJet[nJets] = nCellIn; + nJets++; + } + + } // end of cells loop //reorder jets by et in cone //sort jets by energy Int_t * idx = new Int_t[nJets]; TMath::Sort(nJets, etAlgoJet, idx); - for(Int_t p = 0; p < nJets; p++){ - etaJet[p] = etaAlgoJet[idx[p]]; - phiJet[p] = phiAlgoJet[idx[p]]; - etJet[p] = etAlgoJet[idx[p]]; - etallJet[p] = etAlgoJet[idx[p]]; - ncellsJet[p] = ncellsAlgoJet[idx[p]]; - } - - + for(Int_t p = 0; p < nJets; p++) + { + etaJet[p] = etaAlgoJet[idx[p]]; + phiJet[p] = phiAlgoJet[idx[p]]; + etJet[p] = etAlgoJet[idx[p]]; + etallJet[p] = etAlgoJet[idx[p]]; + ncellsJet[p] = ncellsAlgoJet[idx[p]]; + } + //delete delete index; delete idx; } -//////////////////////////////////////////////////////////////////////// -void AliUA1JetFinderV2::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, - Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, - Int_t* multJet, Int_t* injet) +//////////////////////////////////////////////////////////////////////// +void AliUA1JetFinderV2::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, Float_t* ptT, + Int_t*vectT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* cFlag2T, + Float_t* sFlagT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, + Float_t* etsigJet, Int_t* multJet, Int_t* injet) { - //background subtraction using cone method but without correction in dE/deta distribution - + // + // Background subtraction using cone method but without correction in dE/deta distribution + // Cases to take into account the EMCal geometry are included + // + //calculate energy inside and outside cones AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + Int_t fOpt = fReader->GetReaderHeader()->GetDetector(); Float_t rc= header->GetRadius(); Float_t etIn[30]; Float_t etOut = 0; + + for(Int_t j=0;j<30;j++){etIn[j]=0.;} + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array - // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut - for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; - Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); - if(dr <= rc){ // particles inside this cone - multJet[ijet]++; - injet[jpart] = ijet; - if((fReader->GetCutFlag(jpart)) == 1){ // pt cut - etIn[ijet] += ptT[jpart]; - if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart]; - } - break; - } - }// end jets loop - if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1) - etOut += ptT[jpart]; // particle outside cones and pt cut + + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + multJet[ijet]+=vectT[jpart]; + injet[jpart] = ijet; + + if(cFlagT[jpart] == 1 || cFlag2T[jpart] == 1){ // pt cut + etIn[ijet] += ptT[jpart]; + if(sFlagT[jpart] == 1) etsigJet[ijet]+= ptT[jpart]; + } + break; + } + }// end jets loop + + if(injet[jpart] == -1 && (cFlagT[jpart] == 1 || cFlag2T[jpart] == 1)){ + etOut += ptT[jpart]; // particle outside cones and pt cut + } } //end particle loop //estimate jets and background areas - Float_t areaJet[30]; - Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi(); - for(Int_t k=0; kGetLegoEtaMax())*TMath::Pi(); + + for(Int_t k=0; k header->GetLegoEtaMax()){ // sector outside etamax - Float_t h = header->GetLegoEtaMax() - etaJet[k]; - accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + Float_t h = header->GetLegoEtaMax() - etaJet[k]; + accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); } if(detamin < header->GetLegoEtaMin()){ // sector outside etamin - Float_t h = header->GetLegoEtaMax() + etaJet[k]; - accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + Float_t h = header->GetLegoEtaMax() + etaJet[k]; + accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); } areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin; areaOut = areaOut - areaJet[k]; + } + //subtract background using area method + for(Int_t ljet=0; ljetGetLegoEtaMax())*TMath::Pi(); + etbgTotalN = etOut*areaT/areaOut; + } + else { // If EMCal included + Float_t areaJet[30]; + Float_t areaOut = 2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax() - header->GetLegoPhiMin()); + for(Int_t k=0; kGetLegoEtaMax(); + Float_t eMin = header->GetLegoEtaMin(); + Float_t pMax = header->GetLegoPhiMax(); + Float_t pMin = header->GetLegoPhiMin(); + Float_t accetamax = 0.0; Float_t accetamin = 0.0; + Float_t accphimax = 0.0; Float_t accphimin = 0.0; + if((detamax > eMax && dphimax >= (pMin+2*rc) && dphimax <= pMax )|| + (detamax > eMax && dphimin <= (pMax-2*rc) && dphimin >= pMin )){ + Float_t h = eMax - etaJet[k]; + accetamax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + if((detamin < eMin && dphimax >= (pMin+2*rc) && dphimax <= pMax )|| + (detamin < eMin && dphimin <= (pMax-2*rc) && dphimin >= pMin )){ + Float_t h = eMax + etaJet[k]; + accetamin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + if((dphimax > pMax && detamax >= (eMin+2*rc) && detamax <= eMax )|| + (dphimax > pMax && detamin <= (eMax-2*rc) && detamin >= eMin )){ + Float_t h = pMax - phiJet[k]; + accphimax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + if((dphimin < eMin && detamax >= (eMin+2*rc) && detamax <= eMax )|| + (dphimin < eMin && detamin <= (eMax-2*rc) && detamin >= eMin )){ + Float_t h = phiJet[k] - pMin; + accphimin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + + if(detamax > eMax && dphimax > pMax ){ + Float_t he = eMax - etaJet[k]; + Float_t hp = pMax - phiJet[k]; + Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2)); + Float_t alphae = TMath::ACos(he/rc); + Float_t alphap = TMath::ACos(hp/rc); + Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4; + if(rlim <= rc){ + accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he); + accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp); + } + if(rlim > rc){ + accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he); + accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)- + ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+ + rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad); + } + } + + if(detamax > eMax && dphimin < pMin ){ + Float_t he = eMax - etaJet[k]; + Float_t hp = phiJet[k] - pMin; + Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2)); + Float_t alphae = TMath::ACos(he/rc); + Float_t alphap = TMath::ACos(hp/rc); + Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4; + if(rlim <= rc){ + accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he); + accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp); + } + if(rlim > rc){ + accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he); + accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)- + ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+ + rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad); + } + } + + if(detamin < eMin && dphimax > pMax ){ + Float_t he = eMax + etaJet[k]; + Float_t hp = pMax - phiJet[k]; + Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2)); + Float_t alphae = TMath::ACos(he/rc); + Float_t alphap = TMath::ACos(hp/rc); + Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4; + if(rlim <= rc){ + accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he); + accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp); + } + if(rlim > rc){ + accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he); + accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)- + ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+ + rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad); + } + } + + if(detamin < eMin && dphimin < pMin ){ + Float_t he = eMax + etaJet[k]; + Float_t hp = phiJet[k] - pMin; + Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2)); + Float_t alphae = TMath::ACos(he/rc); + Float_t alphap = TMath::ACos(hp/rc); + Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4; + if(rlim <= rc){ + accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he); + accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp); + } + if(rlim > rc){ + accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he); + accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)- + ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+ + rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad); + } + } + areaJet[k] = rc*rc*TMath::Pi() - accetamax - accetamin - accphimax - accphimin; + areaOut = areaOut - areaJet[k]; + } // end loop on jets + + //subtract background using area method + for(Int_t ljet=0; ljetGetLegoEtaMax()*header->GetLegoPhiMax()); + etbgTotalN = etOut*areaT/areaOut; + } + +} + +//////////////////////////////////////////////////////////////////////// +void AliUA1JetFinderV2::SubtractBackgC(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN, + Float_t* ptT, Float_t* etaT, Float_t* phiT, + Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, + Int_t* multJet, Int_t* injet) +{ + //background subtraction using cone method but without correction in dE/deta distribution + + //calculate energy inside and outside cones + AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; + Float_t rc= header->GetRadius(); + Float_t etIn[30]; + Float_t etOut = 0; + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array + // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + multJet[ijet]++; + injet[jpart] = ijet; + if((fReader->GetCutFlag(jpart)) == 1){ // pt cut + etIn[ijet] += ptT[jpart]; + if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart]; + } + break; + } + }// end jets loop + if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1) + etOut += ptT[jpart]; // particle outside cones and pt cut + } //end particle loop + + //estimate jets and background areas + Float_t areaJet[30]; + Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi(); + for(Int_t k=0; k header->GetLegoEtaMax()){ // sector outside etamax + Float_t h = header->GetLegoEtaMax() - etaJet[k]; + accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + if(detamin < header->GetLegoEtaMin()){ // sector outside etamin + Float_t h = header->GetLegoEtaMax() + etaJet[k]; + accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + } + areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin; + areaOut = areaOut - areaJet[k]; } //subtract background using area method for(Int_t ljet=0; ljetGetLegoEtaMax())*TMath::Pi(); etbgTotalN = etOut*areaT/areaOut; - - + } -//////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////// void AliUA1JetFinderV2::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, + Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet) { @@ -555,295 +1417,291 @@ void AliUA1JetFinderV2::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotal //background subtraction using statistical method AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; Float_t etbgStat = header->GetBackgStat(); // pre-calculated background - + //calculate energy inside Float_t rc= header->GetRadius(); Float_t etIn[30]; - - for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array - //if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut - for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; - Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); - if(dr <= rc){ // particles inside this cone - multJet[ijet]++; - injet[jpart] = ijet; - if((fReader->GetCutFlag(jpart)) == 1){ // pt cut - etIn[ijet]+= ptT[jpart]; - if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; - } - break; - } - }// end jets loop - } //end particle loop - + + for(Int_t jpart = 0; jpart < nIn; jpart++) + { // loop for all particles in array + + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + multJet[ijet]++; + injet[jpart] = ijet; + if(cFlagT[jpart] == 1){ // pt cut + etIn[ijet]+= ptT[jpart]; + if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart]; + } + break; + } + }// end jets loop + } //end particle loop + //calc jets areas Float_t areaJet[30]; Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi(); - for(Int_t k=0; k header->GetLegoEtaMax()){ // sector outside etamax - Float_t h = header->GetLegoEtaMax() - etaJet[k]; - accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + Float_t h = header->GetLegoEtaMax() - etaJet[k]; + accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); } if(detamin < header->GetLegoEtaMin()){ // sector outside etamin - Float_t h = header->GetLegoEtaMax() + etaJet[k]; - accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); + Float_t h = header->GetLegoEtaMax() + etaJet[k]; + accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h); } areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin; - } + } //subtract background using area method for(Int_t ljet=0; ljetGetRadius(); - Float_t etamax = header->GetLegoEtaMax(); - Float_t etamin = header->GetLegoEtaMin(); - Int_t ndiv = 100; - - // jet energy and area arrays - TH1F* hEtJet[30]; - TH1F* hAreaJet[30]; - for(Int_t mjet=0; mjetGetRadius(); + Float_t etamax = header->GetLegoEtaMax(); + Float_t etamin = header->GetLegoEtaMin(); + Int_t ndiv = 100; + + // jet energy and area arrays + TH1F* hEtJet[30]; + TH1F* hAreaJet[30]; + for(Int_t mjet=0; mjet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; - Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); - if(dr <= rc){ // particles inside this cone - injet[jpart] = ijet; - multJet[ijet]++; - if((fReader->GetCutFlag(jpart)) == 1){// pt cut - hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone - if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; - } - break; - } - }// end jets loop - if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1) - hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones + // background energy and area + TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); + TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); + + //fill energies + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + injet[jpart] = ijet; + multJet[ijet]++; + if(cFlagT[jpart] == 1){// pt cut + hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone + if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart]; + } + break; + } + }// end jets loop + + if(injet[jpart] == -1 && cFlagT[jpart] == 1) + hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones } //end particle loop - //calc areas - Float_t eta0 = etamin; - Float_t etaw = (etamax - etamin)/((Float_t)ndiv); - Float_t eta1 = eta0 + etaw; - for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins - Float_t etac = eta0 + etaw/2.0; - Float_t areabg = etaw*2.0*TMath::Pi(); - for(Int_t ijet=0; ijet rc && deta1 < rc){ - acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); - areaj = acc1; - } - if(deta0 < rc && deta1 > rc){ - acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); - areaj = acc0; - } - if(deta0 < rc && deta1 < rc){ - acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); - acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); - if(eta1Fill(etac,areaj); - areabg = areabg - areaj; - } // end jets loop - hAreaBackg->Fill(etac,areabg); - eta0 = eta1; - eta1 = eta1 + etaw; - } // end loop for all eta bins - - //subtract background - for(Int_t kjet=0; kjetGetBinContent(bin)){ - Float_t areab = hAreaBackg->GetBinContent(bin); - Float_t etb = hEtBackg->GetBinContent(bin); - Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab; - etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction - } - } - } - - // calc background total - Double_t etOut = hEtBackg->Integral(); - Double_t areaOut = hAreaBackg->Integral(); - Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi(); - etbgTotalN = etOut*areaT/areaOut; - - //delete - for(Int_t ljet=0; ljet rc && deta1 < rc){ + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + areaj = acc1; + } + if(deta0 < rc && deta1 > rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + areaj = acc0; + } + if(deta0 < rc && deta1 < rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + if(eta1Fill(etac,areaj); + areabg = areabg - areaj; + } // end jets loop + hAreaBackg->Fill(etac,areabg); + eta0 = eta1; + eta1 = eta1 + etaw; + } // end loop for all eta bins + + //subtract background + for(Int_t kjet=0; kjetGetBinContent(bin)){ + Float_t areab = hAreaBackg->GetBinContent(bin); + Float_t etb = hEtBackg->GetBinContent(bin); + Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab; + etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction + } + } + } -//////////////////////////////////////////////////////////////////////// + // calc background total + Double_t etOut = hEtBackg->Integral(); + Double_t areaOut = hAreaBackg->Integral(); + Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi(); + etbgTotalN = etOut*areaT/areaOut; + + //delete + for(Int_t ljet=0; ljetGetBackgCutRatio(); - - - //general - Float_t rc= header->GetRadius(); - Float_t etamax = header->GetLegoEtaMax(); - Float_t etamin = header->GetLegoEtaMin(); - Int_t ndiv = 100; - - // jet energy and area arrays - TH1F* hEtJet[30]; - TH1F* hAreaJet[30]; - for(Int_t mjet=0; mjetGetBackgCutRatio(); + + //general + Float_t rc= header->GetRadius(); + Float_t etamax = header->GetLegoEtaMax(); + Float_t etamin = header->GetLegoEtaMin(); + Int_t ndiv = 100; + + // jet energy and area arrays + TH1F* hEtJet[30]; + TH1F* hAreaJet[30]; + for(Int_t mjet=0; mjetGetCutFlag(jpart)) != 1) continue; - for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; - Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); - if(dr <= rc){ // particles inside this cone - multJet[ijet]++; - injet[jpart] = ijet; - if((fReader->GetCutFlag(jpart)) == 1){ //pt cut - hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut - if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart]; - } - break; - } - }// end jets loop - if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones + // background energy and area + TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax); // change range + TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax); // change range + + //fill energies + for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array + for(Int_t ijet=0; ijet TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi; + Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi); + if(dr <= rc){ // particles inside this cone + multJet[ijet]++; + injet[jpart] = ijet; + if(cFlagT[jpart] == 1){ //pt cut + hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut + if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart]; + } + break; + } + }// end jets loop + if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones } //end particle loop - //calc areas - Float_t eta0 = etamin; - Float_t etaw = (etamax - etamin)/((Float_t)ndiv); - Float_t eta1 = eta0 + etaw; - for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins - Float_t etac = eta0 + etaw/2.0; - Float_t areabg = etaw*2.0*TMath::Pi(); - for(Int_t ijet=0; ijet rc && deta1 < rc){ - acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); - areaj = acc1; - } - if(deta0 < rc && deta1 > rc){ - acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); - areaj = acc0; - } - if(deta0 < rc && deta1 < rc){ - acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); - acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); - if(eta1Fill(etac,areaj); - areabg = areabg - areaj; - } // end jets loop - hAreaBackg->Fill(etac,areabg); - eta0 = eta1; - eta1 = eta1 + etaw; - } // end loop for all eta bins - - //subtract background - for(Int_t kjet=0; kjetGetBinContent(bin)){ - Float_t areab = hAreaBackg->GetBinContent(bin); - Float_t etb = hEtBackg->GetBinContent(bin); - Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab; - etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction - } - } - } - - // calc background total - Double_t etOut = hEtBackg->Integral(); - Double_t areaOut = hAreaBackg->Integral(); - Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi(); - etbgTotalN = etOut*areaT/areaOut; - - //delete - for(Int_t ljet=0; ljet rc && deta1 < rc){ + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + areaj = acc1; + } + if(deta0 < rc && deta1 > rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + areaj = acc0; + } + if(deta0 < rc && deta1 < rc){ + acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0); + acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1); + if(eta1Fill(etac,areaj); + areabg = areabg - areaj; + } // end jets loop + hAreaBackg->Fill(etac,areabg); + eta0 = eta1; + eta1 = eta1 + etaw; + } // end loop for all eta bins + + //subtract background + for(Int_t kjet=0; kjetGetBinContent(bin)){ + Float_t areab = hAreaBackg->GetBinContent(bin); + Float_t etb = hEtBackg->GetBinContent(bin); + Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab; + etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction + } + } + } + + // calc background total + Double_t etOut = hEtBackg->Integral(); + Double_t areaOut = hAreaBackg->Integral(); + Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi(); + etbgTotalN = etOut*areaT/areaOut; + + //delete + for(Int_t ljet=0; ljetReset(); @@ -852,7 +1710,6 @@ void AliUA1JetFinderV2::Reset() } //////////////////////////////////////////////////////////////////////// - void AliUA1JetFinderV2::WriteJHeaderToFile() { AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; @@ -860,22 +1717,22 @@ void AliUA1JetFinderV2::WriteJHeaderToFile() } //////////////////////////////////////////////////////////////////////// - -void AliUA1JetFinderV2::Init() +void AliUA1JetFinderV2::InitTask(TChain* tree) { + // initializes some variables AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader; - // book lego - fLego = new - TH2F("legoH","eta-phi", - header->GetLegoNbinEta(), header->GetLegoEtaMin(), - header->GetLegoEtaMax(), header->GetLegoNbinPhi(), - header->GetLegoPhiMin(), header->GetLegoPhiMax()); - + // book lego + fLego = new TH2F("legoH","eta-phi", + header->GetLegoNbinEta(), header->GetLegoEtaMin(), + header->GetLegoEtaMax(), header->GetLegoNbinPhi(), + header->GetLegoPhiMin(), header->GetLegoPhiMax()); + fDebug = fReader->GetReaderHeader()->GetDebug(); fOpt = fReader->GetReaderHeader()->GetDetector(); - + + // Tasks initialization if(fOpt>0) - fReader->CreateTasks(); + fReader->CreateTasks(tree); } diff --git a/JETAN/AliUA1JetFinderV2.h b/JETAN/AliUA1JetFinderV2.h index b3823458681..eef1429af29 100644 --- a/JETAN/AliUA1JetFinderV2.h +++ b/JETAN/AliUA1JetFinderV2.h @@ -10,12 +10,14 @@ // manages the search for jets // Author: Rafael.Diaz.Valdes@cern.ch // (version in c++) -// Modified to include neutral particles (magali.estienne@ires.in2p3.fr) //--------------------------------------------------------------------- +#include + #include "AliJetFinder.h" class AliUA1JetHeaderV1; class TH2F; +class TChain; class AliUA1JetFinderV2 : public AliJetFinder { @@ -25,35 +27,44 @@ class AliUA1JetFinderV2 : public AliJetFinder ~AliUA1JetFinderV2(); // others + void FindJetsC(); void FindJets(); + void RunAlgoritmC(Float_t EtbgTotal, Double_t dEtTotal, Int_t& nJets, + Float_t* etJet,Float_t* etaJet, Float_t* phiJet, + Float_t* etallJet, Int_t* ncellsJet); void RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell, - Int_t* flagCell, Float_t etbgTotal, Double_t dEtTotal, + Int_t* flagCell, Float_t* etCell2, Float_t* etaCell2, Float_t* phiCell2, + Int_t* flagCell2, Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etallJet, Int_t* ncellsJet); - - - void SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&EtbgTotalN, + + void SubtractBackgC(Int_t& nIn, Int_t&nJ, Float_t&EtbgTotalN, Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,Int_t* multJet, Int_t* injet); + void SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&EtbgTotalN, Float_t* ptT, Int_t* vectT, + Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* cFlag2T, + Float_t* sFlagT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, + Float_t* etsigJet, Int_t* multJet, Int_t* injet); + void SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& EtbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, + Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet); void SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& EtbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, + Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet); void SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&EtbgTotalN, - Float_t* ptT, Float_t* etaT, Float_t* phiT, + Float_t* ptT, Float_t* etaT, Float_t* phiT, Float_t* cFlagT, Float_t* sFlagT, Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet); void Reset(); - void Init(); + void InitTask(TChain* tree); void WriteJHeaderToFile(); protected: -- 2.43.0