From: kleinb Date: Wed, 21 Jan 2009 09:10:28 +0000 (+0000) Subject: Latest changes ny Aria, added CDF charge jet definition, option tp split pos and neg X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=6ef5bfa4abeb538860310c2e65ea1cacd57a7898;p=u%2Fmrichter%2FAliRoot.git Latest changes ny Aria, added CDF charge jet definition, option tp split pos and neg --- diff --git a/PWG4/JetTasks/AliAnalysisTaskUE.cxx b/PWG4/JetTasks/AliAnalysisTaskUE.cxx index 7ca9cc504e7..e00296166aa 100644 --- a/PWG4/JetTasks/AliAnalysisTaskUE.cxx +++ b/PWG4/JetTasks/AliAnalysisTaskUE.cxx @@ -14,7 +14,7 @@ **************************************************************************/ /* $Id:$ */ - + #include #include #include @@ -27,6 +27,8 @@ #include #include #include +#include +#include #include "AliAnalysisTaskUE.h" #include "AliAnalysisManager.h" @@ -68,44 +70,52 @@ ClassImp( AliAnalysisTaskUE) //____________________________________________________________________ AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name): - AliAnalysisTask(name, ""), - fDebug(kFALSE), - fAOD(0x0), - fListOfHistos(0x0), - fBinsPtInHist(30), - fMinJetPtInHist(0.), - fMaxJetPtInHist(300.), - fAnaType(1), - fRegionType(1), - fConeRadius(0.7), - fOrdering(1), - fJet1EtaCut(0.2), - fJet2DeltaPhiCut(2.616), // 150 degrees - fJet2RatioPtCut(0.8), - fJet3PtCut(15.), - fTrackPtCut(0.), - fTrackEtaCut(0.9), - fhNJets(0x0), - fhEleadingPt(0x0), - fhMinRegPtDist(0x0), - fhRegionMultMin(0x0), - fhMinRegAvePt(0x0), - fhMinRegSumPt(0x0), - fhMinRegMaxPtPart(0x0), - fhMinRegSumPtvsMult(0x0), - fhdNdEta_PhiDist(0x0), - fhFullRegPartPtDistVsEt(0x0), - fhTransRegPartPtDistVsEt(0x0), - fhRegionSumPtMaxVsEt(0x0), - fhRegionMultMax(0x0), - fhRegionMultMaxVsEt(0x0), - fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0), - fhRegionMultMinVsEt(0x0), - fhRegionAveSumPtVsEt(0x0), - fhRegionDiffSumPtVsEt(0x0), - fhRegionAvePartPtMaxVsEt(0x0), - fhRegionAvePartPtMinVsEt(0x0), - fhRegionMaxPartPtMaxVsEt(0x0)//, fhValidRegion(0x0) +AliAnalysisTask(name, ""), +fDebug(kFALSE), +fAOD(0x0), +fAODjets(0x0), +fListOfHistos(0x0), +fBinsPtInHist(30), +fMinJetPtInHist(0.), +fMaxJetPtInHist(300.), +fAnaType(1), +fRegionType(1), +fConeRadius(0.7), +fAreaReg(1.5393), // Pi*0.7*0.7 +fUseChPartJet(kFALSE), +fUseSingleCharge(kFALSE), +fUsePositiveCharge(kTRUE), +fOrdering(1), +fFilterBit(0xFF), +fJetsOnFly(kFALSE), +fJet1EtaCut(0.2), +fJet2DeltaPhiCut(2.616), // 150 degrees +fJet2RatioPtCut(0.8), +fJet3PtCut(15.), +fTrackPtCut(0.), +fTrackEtaCut(0.9), +fhNJets(0x0), +fhEleadingPt(0x0), +fhMinRegPtDist(0x0), +fhRegionMultMin(0x0), +fhMinRegAvePt(0x0), +fhMinRegSumPt(0x0), +fhMinRegMaxPtPart(0x0), +fhMinRegSumPtvsMult(0x0), +fhdNdEta_PhiDist(0x0), +fhFullRegPartPtDistVsEt(0x0), +fhTransRegPartPtDistVsEt(0x0), +fhRegionSumPtMaxVsEt(0x0), +fhRegionMultMax(0x0), +fhRegionMultMaxVsEt(0x0), +fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0), +fhRegionMultMinVsEt(0x0), +fhRegionAveSumPtVsEt(0x0), +fhRegionDiffSumPtVsEt(0x0), +fhRegionAvePartPtMaxVsEt(0x0), +fhRegionAvePartPtMinVsEt(0x0), +fhRegionMaxPartPtMaxVsEt(0x0), +fSettingsTree(0x0)//, fhValidRegion(0x0) { // Default constructor // Define input and output slots here @@ -115,218 +125,75 @@ AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name): DefineOutput(0, TList::Class()); } - //____________________________________________________________________ void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/) { -// Connect the input data - -// We need AOD with tracks and jets. -// Since AOD can be either connected to the InputEventHandler or to the OutputEventHandler -// we need to check where it is and get the pointer to AODEvent in the right way + // Connect the input data - if (fDebug > 1) AliInfo("ConnectInputData() \n"); + // We need AOD with tracks and jets. + // Since AOD can be either connected to the InputEventHandler (input chain fron AOD files) + // or to the OutputEventHandler ( AOD is created by a previus task in the train ) + // we need to check where it is and get the pointer to AODEvent in the right way + if (fDebug > 1) AliInfo("ConnectInputData() \n"); TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); if( handler && handler->InheritsFrom("AliAODInputHandler") ) { - fAOD = ((AliAODInputHandler*)handler)->GetEvent(); - } else { + fAOD = ((AliAODInputHandler*)handler)->GetEvent(); + if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler"); + // Case when jets are reconstructed on the fly from AOD tracks + // (the Jet Finder is using the AliJetAODReader) of InputEventHandler + // and put in the OutputEventHandler AOD. Useful whe you want to reconstruct jets with + // different parameters to default ones stored in the AOD or to use a different algorithm + if( fJetsOnFly ) { handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler(); if( handler && handler->InheritsFrom("AliAODHandler") ) { - fAOD = ((AliAODHandler*)handler)->GetAOD(); - } else { - AliFatal("I can't get any AOD Event Handler"); + fAODjets = ((AliAODHandler*)handler)->GetAOD(); + if (fDebug > 1) AliInfo(" ==== Jets from AliAODHandler"); } - } + } else { + fAODjets = fAOD; + if (fDebug > 1) AliInfo(" ==== Jets from AliAODInputHandler"); + } + } else { + handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler(); + if( handler && handler->InheritsFrom("AliAODHandler") ) { + fAOD = ((AliAODHandler*)handler)->GetAOD(); + fAODjets = fAOD; + if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler"); + } else { + AliFatal("I can't get any AOD Event Handler"); + return; + } + } } //____________________________________________________________________ void AliAnalysisTaskUE::CreateOutputObjects() { -// Create the output container -// + // Create the output container + // if (fDebug > 1) AliInfo("CreateOutPutData()"); // // Histograms CreateHistos(); fListOfHistos->SetOwner(kTRUE); -// OpenFile(0); + // OpenFile(0); } //____________________________________________________________________ void AliAnalysisTaskUE::Exec(Option_t */*option*/) { -// Execute analysis for current event -// - AnalyseUE(); + // Execute analysis for current event + // + if ( fDebug > 3 ) AliInfo( " Processing event..." ); + AnalyseUE(); + // Post the data PostData(0, fListOfHistos); } - - -//____________________________________________________________________ -void AliAnalysisTaskUE::Terminate(Option_t */*option*/) -{ -// Terminate analysis -// - - // Normalize histos to region area TODO: - Double_t areafactor = 1.; - if( fRegionType == 1 ) { - areafactor = 1./ (fTrackEtaCut *2. * TMath::Pi()*2.); - } else { - areafactor = 1./ ( fConeRadius * fConeRadius * TMath::Pi() ); - } - - // Draw some Test plot to the screen - if (!gROOT->IsBatch()) { - - // Update pointers reading them from the output slot - fListOfHistos = dynamic_cast (GetOutputData(0)); - if( !fListOfHistos ) { - AliError("Histogram List is not available"); - return; - } - fhEleadingPt = (TH1F*)fListOfHistos->At(1); - fhdNdEta_PhiDist = (TH1F*)fListOfHistos->At(8); - fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11); - fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12); - fhRegionMultMaxVsEt = (TH1F*)fListOfHistos->At(14); - fhRegionMultMinVsEt = (TH1F*)fListOfHistos->At(15); - fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16); - - fhNJets = (TH1F*)fListOfHistos->At(0); - - // Canvas - TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1200,800); - c1->Divide(2,2); - c1->cd(1); - TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); - h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1); - h1r->Scale( areafactor ); - h1r->SetMarkerStyle(20); - h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); - h1r->SetYTitle("P_{T}^{90, max}"); - h1r->DrawCopy("p"); - - c1->cd(2); - - TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); - h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1); - h2r->Scale( areafactor ); - h2r->SetMarkerStyle(20); - h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); - h2r->SetYTitle("P_{T}^{90, min}"); - h2r->DrawCopy("p"); - - c1->cd(3); - TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); - h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1); - h4r->Scale(2.); // make average - h4r->Scale( areafactor ); - h4r->SetYTitle("#DeltaP_{T}^{90}"); - h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); - h4r->SetMarkerStyle(20); - h4r->DrawCopy("p"); - - c1->cd(4); - TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); - TH1F *h6r = new TH1F("hRegionMultMinVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); - - h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1); - h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1); - h5r->Scale( areafactor ); - h6r->Scale( areafactor ); - h5r->SetYTitle("N_{Tracks}^{90}"); - h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); - h5r->SetMarkerStyle(20); - h5r->DrawCopy("p"); - h6r->SetMarkerStyle(21); - h6r->SetMarkerColor(2); - h6r->SetYTitle("N_{Tracks}^{90}"); - h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); - h6r->DrawCopy("p same"); - c1->Update(); - - TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800); - c2->Divide(2,2); - c2->cd(1); - fhEleadingPt->SetMarkerStyle(20); - fhEleadingPt->SetMarkerColor(2); - fhEleadingPt->DrawCopy("p"); - gPad->SetLogy(); - - c2->cd(2); - fhdNdEta_PhiDist->SetMarkerStyle(20); - fhdNdEta_PhiDist->SetMarkerColor(2); - fhdNdEta_PhiDist->DrawCopy("p"); - gPad->SetLogy(); - - c2->cd(3); - fhNJets->DrawCopy(); - - //fhTransRegPartPtDist = (TH1F*)fListOfHistos->At(2); - fhRegionMultMin = (TH1F*)fListOfHistos->At(3); - fhMinRegAvePt = (TH1F*)fListOfHistos->At(4); - fhMinRegSumPt = (TH1F*)fListOfHistos->At(5); - //fhMinRegMaxPtPart = (TH1F*)fListOfHistos->At(6); - fhMinRegSumPtvsMult = (TH1F*)fListOfHistos->At(7); - - // Canvas - TCanvas* c3 = new TCanvas("c3"," p_{T} dist",160,160,1200,800); - c3->Divide(2,2); - c3->cd(1); - /*fhTransRegPartPtDist->SetMarkerStyle(20); - fhTransRegPartPtDist->SetMarkerColor(2); - fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries()); - fhTransRegPartPtDist->DrawCopy("p"); - gPad->SetLogy(); - */ - c3->cd(2); - fhMinRegSumPt->SetMarkerStyle(20); - fhMinRegSumPt->SetMarkerColor(2); - fhMinRegSumPt->Scale(areafactor); - fhMinRegSumPt->DrawCopy("p"); - gPad->SetLogy(); - - c3->cd(3); - fhMinRegAvePt->SetMarkerStyle(20); - fhMinRegAvePt->SetMarkerColor(2); - fhMinRegAvePt->Scale(areafactor); - fhMinRegAvePt->DrawCopy("p"); - gPad->SetLogy(); - - c3->cd(4); - - TH1F *h7r = new TH1F("hRegionMultMinVsMult", "", 21, -0.5, 20.5); - - h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1); - - h7r->SetMarkerStyle(20); - h7r->SetMarkerColor(2); - h7r->DrawCopy("p"); - - c3->Update(); - - -/* c2->cd(3); - fhValidRegion->SetMarkerStyle(20); - fhValidRegion->SetMarkerColor(2); - fhValidRegion->DrawCopy("p"); -*/ - c2->Update(); - } else { - AliInfo(" Batch mode, not histograms will shown..."); - } - - if( fDebug > 1 ) - AliInfo("End analysis"); - -} - //____________________________________________________________________ void AliAnalysisTaskUE::AnalyseUE() @@ -335,6 +202,7 @@ void AliAnalysisTaskUE::AnalyseUE() // Look for correlations on the tranverse regions to // the leading charged jet // + // ------------------------------------------------ // Find Leading Jets 1,2,3 // (could be skipped if Jets are sort by Pt...) @@ -344,82 +212,124 @@ void AliAnalysisTaskUE::AnalyseUE() Int_t index2 = -1; Double_t maxPtJet3 = 0.; // jet 3 need for back to back exclusive Int_t index3 = -1; - Int_t nJets = fAOD->GetNJets(); - for( Int_t i=0; iGetJet(i); - if( jet->Pt() > maxPtJet1 ) { - maxPtJet3 = maxPtJet2; index3 = index2; - maxPtJet2 = maxPtJet1; index2 = index1; - maxPtJet1 = jet->Pt(); index1 = i; - } else if( jet->Pt() > maxPtJet2 ) { - maxPtJet3 = maxPtJet2; index3 = index2; - maxPtJet2 = jet->Pt(); index2 = i; - } else if( jet->Pt() > maxPtJet3 ) { - maxPtJet3 = jet->Pt(); index3 = i; + TVector3 jetVect[3]; + Int_t nJets = 0; + + if( !fUseChPartJet ) { + + // Use AOD Jets + nJets = fAODjets->GetNJets(); + // printf("AOD %d jets \n", nJets); + for( Int_t i=0; iGetJet(i); + Double_t jetPt = jet->Pt();//*1.666; // FIXME Jet Pt Correction ?????!!! + if( jetPt > maxPtJet1 ) { + maxPtJet3 = maxPtJet2; index3 = index2; + maxPtJet2 = maxPtJet1; index2 = index1; + maxPtJet1 = jetPt; index1 = i; + } else if( jetPt > maxPtJet2 ) { + maxPtJet3 = maxPtJet2; index3 = index2; + maxPtJet2 = jetPt; index2 = i; + } else if( jetPt > maxPtJet3 ) { + maxPtJet3 = jetPt; index3 = i; + } + } + if( index1 != -1 ) { + AliAODJet* jet = fAODjets->GetJet(index1); + jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); + } + if( index2 != -1 ) { + AliAODJet* jet = fAODjets->GetJet(index2); + jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); + } + if( index3 != -1 ) { + AliAODJet* jet = fAODjets->GetJet(index3); + jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); + } + + } else { + + // Use "Charged Particle Jets" + TObjArray* jets = FindChargedParticleJets(); + if( jets ) { + nJets = jets->GetEntriesFast(); + if( nJets > 0 ) { + index1 = 0; + AliAODJet* jet = (AliAODJet*)jets->At(0); + maxPtJet1 = jet->Pt(); + jetVect[0].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); + } + if( nJets > 1 ) { + index2 = 1; + AliAODJet* jet = (AliAODJet*)jets->At(1); + maxPtJet1 = jet->Pt(); + jetVect[1].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); + } + if( nJets > 2 ) { + index3 = 2; + AliAODJet* jet = (AliAODJet*)jets->At(2); + maxPtJet1 = jet->Pt(); + jetVect[2].SetXYZ(jet->Px(), jet->Py(), jet->Pz()); + } + + jets->Delete(); + delete jets; } } - - if( index1 < 0 ) { - AliInfo("\nSkipping Event, not jet found..."); - return; - } else - AliInfo(Form("\nPt Leading Jet = %6.1f eta=%5.1f ", maxPtJet1, fAOD->GetJet(index1)->Eta() )); + fhNJets->Fill(nJets); - TVector3 jetVect[2]; + if( fDebug > 1 ) { + if( index1 < 0 ) { + AliInfo("\n Skipping Event, not jet found..."); + return; + } else + AliInfo(Form("\n Pt Leading Jet = %6.1f eta=%5.3f ", maxPtJet1, jetVect[0].Eta() )); + } + // ---------------------------------------------- - // Cut events by jets topology + // Cut events by jets topology // fAnaType: - // 1 = inclusive, + // 1 = inclusive, // - Jet1 |eta| < fJet1EtaCut // 2 = back to back inclusive // - fulfill case 1 // - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut // - Jet2.Pt/Jet1Pt > fJet2RatioPtCut - // 3 = back to back exclusive + // 3 = back to back exclusive // - fulfill case 2 // - Jet3.Pt < fJet3PtCut - Bool_t isInclusive = kFALSE; - - if( TMath::Abs(fAOD->GetJet(index1)->Eta()) > fJet1EtaCut) { - if( fDebug > 1 ) AliInfo("Skipping Event...Jet1 |eta| > fJet1EtaCut"); + if( index1 < 0 || TMath::Abs(jetVect[0].Eta()) > fJet1EtaCut) { + if( fDebug > 1 ) AliInfo("\n Skipping Event...Jet1 |eta| > fJet1EtaCut"); return; } - isInclusive = kTRUE; - jetVect[0].SetXYZ(fAOD->GetJet(index1)->Px(), - fAOD->GetJet(index1)->Py(), - fAOD->GetJet(index1)->Pz()); // back to back inclusive - Bool_t isB2Binclusive = kFALSE; - if( fAnaType > 1 && index2 > 0 && isInclusive) { - jetVect[1].SetXYZ(fAOD->GetJet(index2)->Px(), - fAOD->GetJet(index2)->Py(), - fAOD->GetJet(index2)->Pz()); + if( fAnaType > 1 && index2 == -1 ) { + if( fDebug > 1 ) AliInfo("\n Skipping Event... no second Jet found"); + return; + } + if( fAnaType > 1 && index2 > -1 ) { if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut || - maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) { - if( fDebug > 1 ) AliInfo("Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut"); - return; + maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) { + if( fDebug > 1 ) AliInfo("\n Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut"); + return; } - isB2Binclusive = kTRUE; } - if (isInclusive && !isB2Binclusive && fAnaType>1) return; // back to back exclusive - Bool_t isB2Bexclusive = kFALSE; - if( fAnaType > 2 && index3 > 0 && isB2Binclusive) { - if( fAOD->GetJet(index3)->Pt() > fJet3PtCut ) { - if( fDebug > 1 ) AliInfo("Skipping Event... Jet3.Pt > fJet3PtCut "); + if( fAnaType > 2 && index3 > -1 ) { + if( maxPtJet3 > fJet3PtCut ) { + if( fDebug > 1 ) AliInfo("\n Skipping Event... Jet3.Pt > fJet3PtCut "); return; } - isB2Bexclusive = kTRUE; } - if (isB2Binclusive && !isB2Bexclusive && fAnaType>2) return; - AliInfo(Form("njets = %d",nJets)); fhEleadingPt->Fill( maxPtJet1 ); - + //Area for Normalization Purpose at Display histos + SetRegionArea(jetVect); + // ---------------------------------------------- // Find max and min regions Double_t sumPtRegionPosit = 0.; @@ -432,114 +342,123 @@ void AliAnalysisTaskUE::AnalyseUE() Int_t nTracks = fAOD->GetNTracks(); for (Int_t ipart=0; ipartGetTrack( ipart ); - - if ( !part->Charge() ) continue; + if ( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection + // PID Selection: Reject everything but hadrons + Bool_t isHadron = part->GetMostProbablePID()==AliAODTrack::kPion || + part->GetMostProbablePID()==AliAODTrack::kKaon || + part->GetMostProbablePID()==AliAODTrack::kProton; + if ( !isHadron ) continue; + + if ( !part->Charge() ) continue; //Only charged + if ( fUseSingleCharge ) { // Charge selection + if ( fUsePositiveCharge && part->Charge() < 0.) continue; // keep Positives + if ( !fUsePositiveCharge && part->Charge() > 0.) continue; // keep Negatives + } + if ( part->Pt() < fTrackPtCut ) continue; - + TVector3 partVect(part->Px(), part->Py(), part->Pz()); Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad; if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi(); fhdNdEta_PhiDist->Fill( deltaPhi ); fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); - + Int_t region = IsTrackInsideRegion( jetVect, &partVect ); if (region > 0) { - if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt(); - sumPtRegionPosit += part->Pt(); - nTrackRegionPosit++; - fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); + if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt(); + sumPtRegionPosit += part->Pt(); + nTrackRegionPosit++; + fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); } if (region < 0) { - if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt(); - sumPtRegionNegat += part->Pt(); - nTrackRegionNegat++; - fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); + if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt(); + sumPtRegionNegat += part->Pt(); + nTrackRegionNegat++; + fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); } } - + //How quantities will be sorted before Fill Min and Max Histogram // 1=Plots will be CDF-like // 2=Plots will be Marchesini-like - if (fOrdering==1){ - if (sumPtRegionPosit > sumPtRegionNegat) { - FillSumPtRegion( maxPtJet1, sumPtRegionPosit, sumPtRegionNegat ); + if( fOrdering == 1 ) { + if( sumPtRegionPosit > sumPtRegionNegat ) { + FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg ); } else { - FillSumPtRegion( maxPtJet1, sumPtRegionNegat, sumPtRegionPosit ); + FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg ); } if (nTrackRegionPosit > nTrackRegionNegat ) { - FillMultRegion( maxPtJet1, nTrackRegionPosit, nTrackRegionNegat, sumPtRegionNegat ); + FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg ); } else { - FillMultRegion( maxPtJet1, nTrackRegionNegat, nTrackRegionPosit, sumPtRegionPosit ); + FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg ); } - } else if (fOrdering==2) { + } else if( fOrdering == 2 ) { if (sumPtRegionPosit > sumPtRegionNegat) { - FillSumPtRegion( maxPtJet1, sumPtRegionPosit, sumPtRegionNegat ); - FillMultRegion( maxPtJet1, nTrackRegionPosit, nTrackRegionNegat, sumPtRegionNegat ); + FillSumPtRegion( maxPtJet1, sumPtRegionPosit/fAreaReg, sumPtRegionNegat/fAreaReg ); + FillMultRegion( maxPtJet1, nTrackRegionPosit/fAreaReg, nTrackRegionNegat/fAreaReg, sumPtRegionNegat/fAreaReg ); } else { - FillSumPtRegion( maxPtJet1, sumPtRegionNegat, sumPtRegionPosit ); - FillMultRegion( maxPtJet1, nTrackRegionNegat, nTrackRegionPosit, sumPtRegionPosit ); + FillSumPtRegion( maxPtJet1, sumPtRegionNegat/fAreaReg, sumPtRegionPosit/fAreaReg ); + FillMultRegion( maxPtJet1, nTrackRegionNegat/fAreaReg, nTrackRegionPosit/fAreaReg, sumPtRegionPosit/fAreaReg ); } } Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.; Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.; - if (avePosRegion > aveNegRegion) { - FillAvePartPtRegion( maxPtJet1, avePosRegion, aveNegRegion ); + if( avePosRegion > aveNegRegion ) { + FillAvePartPtRegion( maxPtJet1, avePosRegion/fAreaReg, aveNegRegion/fAreaReg ); } else { - FillAvePartPtRegion( maxPtJet1, aveNegRegion, avePosRegion ); + FillAvePartPtRegion( maxPtJet1, aveNegRegion/fAreaReg, avePosRegion/fAreaReg ); } - + fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion ); - - // Compute pedestal like magnitude - fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/2.0); - fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/2.0); - + + // Compute pedestal like magnitudes + fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/(2.0*fAreaReg)); + fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/(2.0*fAreaReg)); + } - //____________________________________________________________________ void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ) { - // Fill sumPt of control regions - - // Max cone - fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax ); - // Min cone - fhRegionSumPtMinVsEt->Fill( leadingE, ptMin ); - // MAke distributions for UE comparison with MB data - fhMinRegSumPt->Fill(ptMin); + // Fill sumPt of control regions + + // Max cone + fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax ); + // Min cone + fhRegionSumPtMinVsEt->Fill( leadingE, ptMin ); + // MAke distributions for UE comparison with MB data + fhMinRegSumPt->Fill(ptMin); } //____________________________________________________________________ void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ) { - // Fill average particle Pt of control regions - - // Max cone - fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax ); - // Min cone - fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin ); - // MAke distributions for UE comparison with MB data - fhMinRegAvePt->Fill(ptMin); + // Fill average particle Pt of control regions + + // Max cone + fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax ); + // Min cone + fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin ); + // MAke distributions for UE comparison with MB data + fhMinRegAvePt->Fill(ptMin); } //____________________________________________________________________ void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin ) { - // Fill Nch multiplicity of control regions - - // Max cone - fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax ); - fhRegionMultMax->Fill( nTrackPtmax ); - // Min cone - fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin ); - fhRegionMultMin->Fill( nTrackPtmin ); - // MAke distributions for UE comparison with MB data - fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin); - + // Fill Nch multiplicity of control regions + + // Max cone + fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax ); + fhRegionMultMax->Fill( nTrackPtmax ); + // Min cone + fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin ); + fhRegionMultMin->Fill( nTrackPtmin ); + // MAke distributions for UE comparison with MB data + fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin); } //____________________________________________________________________ @@ -558,7 +477,7 @@ Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVe // transverse regions if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1; if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1; - + } else if( fRegionType == 2 ) { // Cone regions Double_t deltaR = 0.; @@ -568,24 +487,185 @@ Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVe negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2()); if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) { - region = 1; - deltaR = positVect.DrEtaPhi(*partVect); + region = 1; + deltaR = positVect.DrEtaPhi(*partVect); } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) { - region = -1; - deltaR = negatVect.DrEtaPhi(*partVect); + region = -1; + deltaR = negatVect.DrEtaPhi(*partVect); } if (deltaR > fConeRadius) region = 0; - + } else AliError("Unknow region type"); - + // For debug (to be removed) //if( region != 0 ) fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) ); return region; } + +//____________________________________________________________________ +TObjArray* AliAnalysisTaskUE::FindChargedParticleJets() +{ + // Return a TObjArray of "charged particle jets" + // + // Charged particle jet definition from reference: + // + // "Charged jet evolution and the underlying event + // in proton-antiproton collisions at 1.8 TeV" + // PHYSICAL REVIEW D 65 092002, CDF Collaboration + // + // We define "jets" as circular regions in eta-phi space with + // radius defined by R = sqrt( (eta-eta0)^2 +(phi-phi0)^2 ). + // Our jet algorithm is as follows: + // 1- Order all charged particles according to their pT . + // 2- Start with the highest pT particle and include in the jet all + // particles within the radius R=0.7 considering each particle + // in the order of decreasing pT and recalculating the centroid + // of the jet after each new particle is added to the jet . + // 3- Go to the next highest pT particle not already included in + // a jet and add to the jet all particles not already included in + // a jet within R=0.7. + // 4- Continue until all particles are in a jet. + // We define the transverse momentum of the jet to be + // the scalar pT sum of all the particles within the jet, where pT + // is measured with respect to the beam axis + + // 1 - Order all charged particles according to their pT . + Int_t nTracks = fAOD->GetNTracks(); + if( !nTracks ) return 0; + TObjArray tracks(nTracks); + + for (Int_t ipart=0; ipartGetTrack( ipart ); + if( !part->TestFilterBit(fFilterBit) ) continue; // track cut selection + if( !part->Charge() ) continue; + if( part->Pt() < fTrackPtCut ) continue; + tracks.AddLast(part); + } + QSortTracks( tracks, 0, tracks.GetEntriesFast() ); + + nTracks = tracks.GetEntriesFast(); + if( !nTracks ) return 0; + TObjArray *jets = new TObjArray(nTracks); + TIter itrack(&tracks); + while( nTracks ) { + // 2- Start with the highest pT particle ... + Float_t px,py,pz,pt; + AliAODTrack* track = (AliAODTrack*)itrack.Next(); + if( !track ) continue; + px = track->Px(); + py = track->Py(); + pz = track->Pz(); + pt = track->Pt(); // Use the energy member to store Pt + jets->AddLast( new TLorentzVector(px, py, pz, pt) ); + tracks.Remove( track ); + TLorentzVector* jet = (TLorentzVector*)jets->Last(); + // 3- Go to the next highest pT particle not already included... + AliAODTrack* track1; + while ( (track1 = (AliAODTrack*)(itrack.Next())) ) { + Double_t dphi = TVector2::Phi_mpi_pi(jet->Phi()-track1->Phi()); + Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) + + dphi*dphi ); + if( r < fConeRadius ) { + Double_t Pt = jet->E()+track1->Pt(); // Scalar sum of Pt + // recalculating the centroid + Double_t eta = jet->Eta()*jet->E()/Pt + track1->Eta()/Pt; + Double_t phi = jet->Phi()*jet->E()/Pt + track1->Phi()/Pt; + jet->SetPtEtaPhiE( 1., eta, phi, Pt ); + tracks.Remove( track1 ); + } + } + + tracks.Compress(); + nTracks = tracks.GetEntries(); + // 4- Continue until all particles are in a jet. + itrack.Reset(); + } // end while nTracks + + // Convert to AODjets.... + Int_t njets = jets->GetEntriesFast(); + TObjArray* aodjets = new TObjArray(njets); + aodjets->SetOwner(kTRUE); + for(Int_t ijet=0; ijetAt(ijet); + Float_t px, py,pz,en; // convert to 4-vector + px = jet->E() * TMath::Cos(jet->Phi()); // Pt * cos(phi) + py = jet->E() * TMath::Sin(jet->Phi()); // Pt * sin(phi) + pz = jet->E() / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-jet->Eta()))); + en = TMath::Sqrt(px * px + py * py + pz * pz); + aodjets->AddLast( new AliAODJet(px, py, pz, en) ); + } + // Order jets according to their pT . + QSortTracks( *aodjets, 0, aodjets->GetEntriesFast() ); + + // debug + if (fDebug>3) AliInfo(Form(" %d Charged jets found\n",njets)); + + return aodjets; +} + +//____________________________________________________________________ +void AliAnalysisTaskUE::QSortTracks(TObjArray &a, Int_t first, Int_t last) +{ + // Sort array of TObjArray of tracks by Pt using a quicksort algorithm. + + static TObject *tmp; + static int i; // "static" to save stack space + int j; + + while (last - first > 1) { + i = first; + j = last; + for (;;) { + while (++i < last && ((AliVParticle*)a[i])->Pt() > ((AliVParticle*)a[first])->Pt() ) + ; + while (--j > first && ((AliVParticle*)a[j])->Pt() < ((AliVParticle*)a[first])->Pt() ) + ; + if (i >= j) + break; + + tmp = a[i]; + a[i] = a[j]; + a[j] = tmp; + } + if (j == first) { + ++first; + continue; + } + tmp = a[first]; + a[first] = a[j]; + a[j] = tmp; + if (j - first < last - (j + 1)) { + QSortTracks(a, first, j); + first = j + 1; // QSortTracks(j + 1, last); + } else { + QSortTracks(a, j + 1, last); + last = j; // QSortTracks(first, j); + } + } +} + +//____________________________________________________________________ +void AliAnalysisTaskUE::SetRegionArea(TVector3 *jetVect) +{ + Double_t fAreaCorrFactor=0.; + Double_t deltaEta = 0.; + if (fRegionType==1) fAreaReg = 2.*fTrackEtaCut*60.*TMath::Pi()/180.; + else if (fRegionType==2){ + deltaEta = 0.9-TMath::Abs(jetVect[0].Eta()); + if (deltaEta>fConeRadius) fAreaReg = TMath::Pi()*fConeRadius*fConeRadius; + else{ + fAreaCorrFactor = fConeRadius*fConeRadius*TMath::ACos(deltaEta/fConeRadius) - + (deltaEta*fConeRadius)*TMath::Sqrt( 1. - deltaEta*deltaEta/(fConeRadius*fConeRadius)); + fAreaReg=TMath::Pi()*fConeRadius*fConeRadius-fAreaCorrFactor; + } + }else AliWarning("Unknown Rgion Type"); + if (fDebug>5) AliInfo(Form("\n dEta=%5.3f Angle =%5.3f Region Area = %5.3f Corr Factor=%5.4f \n",deltaEta,TMath::ACos(deltaEta/fConeRadius),fAreaReg,fAreaCorrFactor)); +} + //____________________________________________________________________ void AliAnalysisTaskUE::CreateHistos() { @@ -595,7 +675,6 @@ void AliAnalysisTaskUE::CreateHistos() fhNJets = new TH1F("fhNJets", "n Jet", 10, 0, 10); fhNJets->SetXTitle("# of jets"); fhNJets->Sumw2(); - fhNJets->Sumw2(); fListOfHistos->Add( fhNJets ); // At(0) fhEleadingPt = new TH1F("hEleadingPt", "Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); @@ -625,7 +704,7 @@ void AliAnalysisTaskUE::CreateHistos() fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)"); fhMinRegSumPt->Sumw2(); fListOfHistos->Add( fhMinRegSumPt ); // At(5) - + fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ", 50, 0., 20.); fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})"); fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)"); @@ -637,24 +716,24 @@ void AliAnalysisTaskUE::CreateHistos() fhMinRegSumPtvsMult->SetXTitle("N_{charge}"); fhMinRegSumPtvsMult->Sumw2(); fListOfHistos->Add( fhMinRegSumPtvsMult ); // At(7); - + fhdNdEta_PhiDist = new TH1F("hdNdEta_PhiDist", "Charge particle density |#eta|< 1 vs #Delta#phi", 120, 0., 2.*TMath::Pi()); fhdNdEta_PhiDist->SetXTitle("#Delta#phi"); fhdNdEta_PhiDist->SetYTitle("dN_{ch}/d#etad#phi"); fhdNdEta_PhiDist->Sumw2(); fListOfHistos->Add( fhdNdEta_PhiDist ); // At(8) - + // Can be use to get part pt distribution for differente Jet Pt bins fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", "dN/dP_{T} |#eta|<1 vs Leading Jet P_{T}", - 50,0.,50., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + 50,0.,50., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}"); fhFullRegPartPtDistVsEt->SetXTitle("p_{T}"); fhFullRegPartPtDistVsEt->Sumw2(); fListOfHistos->Add( fhFullRegPartPtDistVsEt ); // At(9) - - // Can be use to get part pt distribution for differente Jet Pt bins + + // Can be use to get part pt distribution for differente Jet Pt bins fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", "dN/dP_{T} in tranvese regions |#eta|<1 vs Leading Jet P_{T}", - 50,0.,50., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + 50,0.,50., fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}"); fhTransRegPartPtDistVsEt->SetXTitle("p_{T}"); fhTransRegPartPtDistVsEt->Sumw2(); @@ -665,7 +744,7 @@ void AliAnalysisTaskUE::CreateHistos() fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)"); fhRegionSumPtMaxVsEt->Sumw2(); fListOfHistos->Add( fhRegionSumPtMaxVsEt ); // At(11) - + fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt", "P_{T}^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)"); fhRegionSumPtMinVsEt->Sumw2(); @@ -675,22 +754,22 @@ void AliAnalysisTaskUE::CreateHistos() fhRegionMultMax->SetXTitle("N_{ch tracks}"); fhRegionMultMax->Sumw2(); fListOfHistos->Add( fhRegionMultMax ); // At(13) - + fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt", "N_{ch}^{90, max} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)"); fhRegionMultMaxVsEt->Sumw2(); fListOfHistos->Add( fhRegionMultMaxVsEt ); // At(14) - + fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt", "N_{ch}^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); fhRegionMultMinVsEt->SetXTitle("E (GeV/c)"); fhRegionMultMinVsEt->Sumw2(); fListOfHistos->Add( fhRegionMultMinVsEt ); // At(15) - + fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)"); fhRegionAveSumPtVsEt->Sumw2(); fListOfHistos->Add( fhRegionAveSumPtVsEt ); // At(16) - + fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)"); fhRegionDiffSumPtVsEt->Sumw2(); @@ -700,26 +779,209 @@ void AliAnalysisTaskUE::CreateHistos() fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)"); fhRegionAvePartPtMaxVsEt->Sumw2(); fListOfHistos->Add( fhRegionAvePartPtMaxVsEt ); // At(18) - + fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)"); fhRegionAvePartPtMinVsEt->Sumw2(); fListOfHistos->Add( fhRegionAvePartPtMinVsEt ); // At(19) - + fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)"); fhRegionMaxPartPtMaxVsEt->Sumw2(); fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt ); // At(20) - -/* - // For debug region selection - fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi", - 20, -2.,2., 62, -TMath::Pi(), TMath::Pi()); - fhValidRegion->SetYTitle("#Delta#phi"); - fhValidRegion->Sumw2(); - fListOfHistos->Add( fhValidRegion ); // At(15) -*/ + + fSettingsTree = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation"); + fListOfHistos->Add( fSettingsTree ); // At(21) + + /* + // For debug region selection + fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi", + 20, -2.,2., 62, -TMath::Pi(), TMath::Pi()); + fhValidRegion->SetYTitle("#Delta#phi"); + fhValidRegion->Sumw2(); + fListOfHistos->Add( fhValidRegion ); // At(15) + */ } - +//____________________________________________________________________ +void AliAnalysisTaskUE::Terminate(Option_t */*option*/) +{ + // Terminate analysis + // + + //Save Analysis Settings + WriteSettings(); + + // Normalize histos to region area TODO: + // Normalization done at Analysis, taking into account + // area variations on per-event basis (cone case) + + // Draw some Test plot to the screen + if (!gROOT->IsBatch()) { + + // Update pointers reading them from the output slot + fListOfHistos = dynamic_cast (GetOutputData(0)); + if( !fListOfHistos ) { + AliError("Histogram List is not available"); + return; + } + fhEleadingPt = (TH1F*)fListOfHistos->At(1); + fhdNdEta_PhiDist = (TH1F*)fListOfHistos->At(8); + fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11); + fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12); + fhRegionMultMaxVsEt = (TH1F*)fListOfHistos->At(14); + fhRegionMultMinVsEt = (TH1F*)fListOfHistos->At(15); + fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16); + + fhNJets = (TH1F*)fListOfHistos->At(0); + + // Canvas + TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1200,800); + c1->Divide(2,2); + c1->cd(1); + TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1); + //h1r->Scale( areafactor ); + h1r->SetMarkerStyle(20); + h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); + h1r->SetYTitle("P_{T}^{90, max}"); + h1r->DrawCopy("p"); + + c1->cd(2); + + TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1); + //h2r->Scale( areafactor ); + h2r->SetMarkerStyle(20); + h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); + h2r->SetYTitle("P_{T}^{90, min}"); + h2r->DrawCopy("p"); + + c1->cd(3); + TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1); + h4r->Scale(2.); // make average + //h4r->Scale( areafactor ); + h4r->SetYTitle("#DeltaP_{T}^{90}"); + h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); + h4r->SetMarkerStyle(20); + h4r->DrawCopy("p"); + + c1->cd(4); + TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + TH1F *h6r = new TH1F("hRegionMultMinVsEtleading", "", fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist); + + h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1); + h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1); + //h5r->Scale( areafactor ); + //h6r->Scale( areafactor ); + h5r->SetYTitle("N_{Tracks}^{90}"); + h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); + h5r->SetMarkerStyle(20); + h5r->DrawCopy("p"); + h6r->SetMarkerStyle(21); + h6r->SetMarkerColor(2); + h6r->SetYTitle("N_{Tracks}^{90}"); + h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)"); + h6r->DrawCopy("p same"); + c1->Update(); + + TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800); + c2->Divide(2,2); + c2->cd(1); + fhEleadingPt->SetMarkerStyle(20); + fhEleadingPt->SetMarkerColor(2); + fhEleadingPt->DrawCopy("p"); + gPad->SetLogy(); + + c2->cd(2); + fhdNdEta_PhiDist->SetMarkerStyle(20); + fhdNdEta_PhiDist->SetMarkerColor(2); + fhdNdEta_PhiDist->DrawCopy("p"); + gPad->SetLogy(); + + c2->cd(3); + fhNJets->DrawCopy(); + + //fhTransRegPartPtDist = (TH1F*)fListOfHistos->At(2); + fhRegionMultMin = (TH1F*)fListOfHistos->At(3); + fhMinRegAvePt = (TH1F*)fListOfHistos->At(4); + fhMinRegSumPt = (TH1F*)fListOfHistos->At(5); + //fhMinRegMaxPtPart = (TH1F*)fListOfHistos->At(6); + fhMinRegSumPtvsMult = (TH1F*)fListOfHistos->At(7); + + // Canvas + TCanvas* c3 = new TCanvas("c3"," p_{T} dist",160,160,1200,800); + c3->Divide(2,2); + c3->cd(1); + /*fhTransRegPartPtDist->SetMarkerStyle(20); + fhTransRegPartPtDist->SetMarkerColor(2); + fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries()); + fhTransRegPartPtDist->DrawCopy("p"); + gPad->SetLogy(); + */ + c3->cd(2); + fhMinRegSumPt->SetMarkerStyle(20); + fhMinRegSumPt->SetMarkerColor(2); + //fhMinRegSumPt->Scale(areafactor); + fhMinRegSumPt->DrawCopy("p"); + gPad->SetLogy(); + + c3->cd(3); + fhMinRegAvePt->SetMarkerStyle(20); + fhMinRegAvePt->SetMarkerColor(2); + //fhMinRegAvePt->Scale(areafactor); + fhMinRegAvePt->DrawCopy("p"); + gPad->SetLogy(); + + c3->cd(4); + + TH1F *h7r = new TH1F("hRegionMultMinVsMult", "", 21, -0.5, 20.5); + + h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1); + + h7r->SetMarkerStyle(20); + h7r->SetMarkerColor(2); + h7r->DrawCopy("p"); + + c3->Update(); + + + /* c2->cd(3); + fhValidRegion->SetMarkerStyle(20); + fhValidRegion->SetMarkerColor(2); + fhValidRegion->DrawCopy("p"); + */ + c2->Update(); + } else { + AliInfo(" Batch mode, not histograms will be shown..."); + } + + if( fDebug > 1 ) + AliInfo("End analysis"); + +} +void AliAnalysisTaskUE::WriteSettings() +{ + + fSettingsTree->Branch("fAnaType", &fAnaType, "Ana/I"); + fSettingsTree->Branch("fRegionType", &fRegionType,"Reg/I"); + fSettingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D"); + fSettingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O"); + fSettingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O"); + fSettingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O"); + fSettingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I"); + fSettingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D"); + fSettingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D"); + fSettingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D"); + fSettingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D"); + fSettingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D"); + fSettingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D"); + fSettingsTree->Fill(); + + if (fDebug>5){ + AliInfo(" All Analysis Settings in Saved Tree"); + fSettingsTree->Scan(); + } +} diff --git a/PWG4/JetTasks/AliAnalysisTaskUE.h b/PWG4/JetTasks/AliAnalysisTaskUE.h index 0e3ad6274c3..5a8a904611e 100644 --- a/PWG4/JetTasks/AliAnalysisTaskUE.h +++ b/PWG4/JetTasks/AliAnalysisTaskUE.h @@ -1,9 +1,9 @@ #ifndef ALIANALYSISTASKUE_H #define ALIANALYSISTASKUE_H - + /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ - + #include "AliAnalysisTaskSE.h" class AliESDEvent; @@ -12,127 +12,158 @@ class TH1F; class TH2F; class TH1I; class TVector3; - -class AliAnalysisTaskUE : public AliAnalysisTask -{ - public: - AliAnalysisTaskUE(const char* name="AliAnalysisTaskUE"); - virtual ~AliAnalysisTaskUE() {;} - - // Implementation of interface methods - virtual void ConnectInputData(Option_t *); - virtual void CreateOutputObjects(); - virtual void Exec(Option_t *option); - virtual void Terminate(Option_t *); - - // Setters - virtual void SetDebugLevel(Int_t level) { fDebug = level; } - void SetPtRangeInHist(Int_t bin, Double_t min, Double_t max) { - fBinsPtInHist = bin; - fMinJetPtInHist = min; - fMaxJetPtInHist = max; - } - - void SetAnaTopology(Int_t val) { fAnaType = val; } - void SetRegionType(Int_t val) { fRegionType = val; } - void SetConeRadius(Double_t val) { fConeRadius = val; } - void SetPtSumOrdering(Int_t val) { fOrdering = val; } - // Jet cuts - void SetJet1EtaCut(Double_t val) { fJet1EtaCut = val; } - void SetJet2DeltaPhiCut(Double_t val) { fJet2DeltaPhiCut = val; } - void SetJet2RatioPtCut(Double_t val) { fJet2RatioPtCut = val; } - void SetJet3PtCut(Double_t val) { fJet3PtCut = val; } - // track cuts - void SetTrackPtCut(Double_t val) { fTrackPtCut = val; } - void SetTrackEtaCut(Double_t val) { fTrackEtaCut = val; } - -private: - AliAnalysisTaskUE(const AliAnalysisTaskUE &det); -AliAnalysisTaskUE& operator=(const AliAnalysisTaskUE &det); +class TTree; - void AnalyseUE(); - Int_t IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect); - void CreateHistos(); - void FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ); - void FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ); - void FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin ); - - - Int_t fDebug; // Debug flag - AliAODEvent* fAOD; //! AOD Event - TList* fListOfHistos; // Output list of histograms - - // Config - Int_t fBinsPtInHist; // # bins for Pt histos range - Double_t fMinJetPtInHist; // min Jet Pt value for histo range - Double_t fMaxJetPtInHist; // max Jet Pt value for histo range - - // Cuts - Int_t fAnaType; // Analysis type on jet topology: - // 1=inclusive (default) - // 2=back to back inclusive - // 3=back to back exclusive - // 4=gama jet (back to back) ??? - // Minimum bias - // 31 = Semi jet (charged leading particle jets) - // 32 = Random jetcone ? - // 33 = Swiss chees ? - - Int_t fRegionType; // 1 = transverse regions (default) - // 2 = cone regions - - Double_t fConeRadius; // if selected Cone-like region type set Radius (=0.7 default) - - Int_t fOrdering; // Pt and multiplicity summation ordering: - // 1=CDF-like -independent sorting according quantity to be scored: Double sorting- (default) - // if Pt summation will be scored take Pt minimum between both zones and - // fill Pt Max. and Min. histog. accordingly - // if Multiplicity summation will be scored take Mult. minimum between both zones and - // fill Mult Max and Min histog. accordingly - // 2=Marchesini-like (Only Pt sorting: Single sorting) - // sort only according Pt summation scored, find minimum between both zones and - // fill Pt and Multiplicity Max and Min summation histog. following only this criterium - // 3=User Selection sorting (NOTE: USER must implement it within cxx) - - // Jet cuts - Double_t fJet1EtaCut; // |jet1 eta| < fJet1EtaCut (fAnaType = 1,2,3) - Double_t fJet2DeltaPhiCut; // |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut (fAnaType = 2,3) - Double_t fJet2RatioPtCut; // Jet2.Pt/Jet1Pt > fJet2RatioPtCut (fAnaType = 2,3) - Double_t fJet3PtCut; // Jet3.Pt < fJet3PtCut (fAnaType = 3) - // track cuts - Double_t fTrackPtCut; // Pt cut of tracks in the regions - Double_t fTrackEtaCut; // Eta cut on tracks in the regions (fRegionType=1) +class AliAnalysisTaskUE : public AliAnalysisTask + { + public: + AliAnalysisTaskUE(const char* name="AliAnalysisTaskUE"); + virtual ~AliAnalysisTaskUE() {;} + + // Implementation of interface methods + virtual void ConnectInputData(Option_t *); + virtual void CreateOutputObjects(); + virtual void Exec(Option_t *option); + virtual void Terminate(Option_t *); + + // Setters + virtual void SetDebugLevel( Int_t level ) { fDebug = level; } + void SetPtRangeInHist( Int_t bin, Double_t min, Double_t max ) { + fBinsPtInHist = bin; + fMinJetPtInHist = min; + fMaxJetPtInHist = max; + } + + void SetAnaTopology( Int_t val ) { fAnaType = val; } + void SetRegionType( Int_t val ) { fRegionType = val; } + void SetUseChPartJet( Int_t val ) { fUseChPartJet = val; } + void SetPtSumOrdering( Bool_t val ) { fOrdering = val; } + void SetFilterBit( UInt_t val ) { fFilterBit = val; } + void SetJetsOnFly( Bool_t val ) { fJetsOnFly = val; } + void SetConeRadius( Double_t val ) { fConeRadius = val; } + void SetUseSingleCharge() { fUseSingleCharge = kTRUE; } + void SetUseNegativeChargeType() { fUsePositiveCharge = kFALSE; } + // Jet cuts + void SetJet1EtaCut( Double_t val ) { fJet1EtaCut = val; } + void SetJet2DeltaPhiCut( Double_t val ) { fJet2DeltaPhiCut = val; } + void SetJet2RatioPtCut( Double_t val ) { fJet2RatioPtCut = val; } + void SetJet3PtCut( Double_t val ) { fJet3PtCut = val; } + // track cuts + void SetTrackPtCut( Double_t val ) { fTrackPtCut = val; } + void SetTrackEtaCut( Double_t val ) { fTrackEtaCut = val; } + + private: + AliAnalysisTaskUE(const AliAnalysisTaskUE &det); + AliAnalysisTaskUE& operator=(const AliAnalysisTaskUE &det); + + void AnalyseUE(); + Int_t IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect); + void CreateHistos(); + void SetRegionArea(TVector3 *jetVect); + void FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ); + void FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin ); + void FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin ); + TObjArray* FindChargedParticleJets(); + void QSortTracks(TObjArray &a, Int_t first, Int_t last); + void WriteSettings(); + + Int_t fDebug; // Debug flag + AliAODEvent* fAOD; //! AOD Event + AliAODEvent* fAODjets; //! AOD Event for reconstructed on the fly (see ConnectInputData() + TList* fListOfHistos; // Output list of histograms + + // Config + Int_t fBinsPtInHist; // # bins for Pt histos range + Double_t fMinJetPtInHist; // min Jet Pt value for histo range + Double_t fMaxJetPtInHist; // max Jet Pt value for histo range + + // Cuts + Int_t fAnaType; // Analysis type on jet topology: + // 1=inclusive (default) + // 2=back to back inclusive + // 3=back to back exclusive + // 4=gama jet (back to back) ??? + // Minimum bias + // 31 = Semi jet (charged leading particle jets) + // 32 = Random jetcone ? + // 33 = Swiss chees ? + + // UE analysis is conducted in different type of regions + // Transverse are those like defined in: R. Field Acta Physica Polonica B. Vol 36 No. 2 pg 167 (2005) + // Cone regions like defined in: Phys. Rev. D 70, 072002 (2004) + Int_t fRegionType; // 1 = transverse regions (default) + // 2 = cone regions + Double_t fConeRadius; // if selected Cone-like region type, set Radius (=0.7 default) + Double_t fAreaReg; // Area of the region To be used as normalization factor when filling histograms + // if fRegionType = 2 not always it is included within eta range + Bool_t fUseChPartJet; // Use "Charged Particle Jet" instead of jets from AOD + // see FindChargedParticleJets() + + // Theoreticians ask for tools charge-aware + // especially those which are related to multiplicity and MC-tunings + // see arXiv:hep-ph/0507008v3 + Bool_t fUseSingleCharge; //Make analysis for a single type of charge (=kFALSE default) + Bool_t fUsePositiveCharge; //If Single type of charge used then set which one (=kTRUE default positive) + + Int_t fOrdering; // Pt and multiplicity summation ordering: + // 1=CDF-like -independent sorting according quantity to be scored: Double sorting- (default) + // if Pt summation will be scored take Pt minimum between both zones and + // fill Pt Max. and Min. histog. accordingly + // if Multiplicity summation will be scored take Mult. minimum between both zones and + // fill Mult Max and Min histog. accordingly + // Bib: + // 2=Marchesini-like (Only Pt sorting: Single sorting) + // sort only according Pt summation scored, find minimum between both zones and + // fill Pt and Multiplicity Max and Min summation histog. following only this criterium + // Bib: Phys. Rev. D 38, 3419 (1988) + // 3=User Selection sorting (NOTE: USER must implement it within cxx) + + UInt_t fFilterBit; // Select tracks from an specific track cut (default 0xFF all track selected) + Bool_t fJetsOnFly; // if jets are reconstructed on the fly from AOD tracks (see ConnectInputData() ) + + // Jet cuts + Double_t fJet1EtaCut; // |jet1 eta| < fJet1EtaCut (fAnaType = 1,2,3) + Double_t fJet2DeltaPhiCut; // |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut (fAnaType = 2,3) + Double_t fJet2RatioPtCut; // Jet2.Pt/Jet1Pt > fJet2RatioPtCut (fAnaType = 2,3) + Double_t fJet3PtCut; // Jet3.Pt < fJet3PtCut (fAnaType = 3) + // track cuts + Double_t fTrackPtCut; // Pt cut of tracks in the regions + Double_t fTrackEtaCut; // Eta cut on tracks in the regions (fRegionType=1) + + // Histograms ( are owned by fListOfHistos TList ) + TH1F* fhNJets; //! + TH1F* fhEleadingPt; //! + + TH1F* fhMinRegPtDist; //! + TH1F* fhRegionMultMin; //! + TH1F* fhMinRegAvePt; //! + TH1F* fhMinRegSumPt; //! + TH1F* fhMinRegMaxPtPart; //! + TH1F* fhMinRegSumPtvsMult; //! + + TH1F* fhdNdEta_PhiDist; //! + TH2F* fhFullRegPartPtDistVsEt; //! + TH2F* fhTransRegPartPtDistVsEt; //! + + TH1F* fhRegionSumPtMaxVsEt; //! + TH1I* fhRegionMultMax; //! + TH1F* fhRegionMultMaxVsEt; //! + TH1F* fhRegionSumPtMinVsEt; //! + TH1F* fhRegionMultMinVsEt; //! + TH1F* fhRegionAveSumPtVsEt; //! + TH1F* fhRegionDiffSumPtVsEt; //! + + TH1F* fhRegionAvePartPtMaxVsEt; //! + TH1F* fhRegionAvePartPtMinVsEt; //! + TH1F* fhRegionMaxPartPtMaxVsEt; //! + + TTree* fSettingsTree; //! Fast Settings saving + + // TH2F* fhValidRegion; //! test to be canceled + + ClassDef( AliAnalysisTaskUE, 1); // Analysis task for Underlying Event analysis + }; - // Histograms ( are owner by fListOfHistos TList ) - TH1F* fhNJets; //! - TH1F* fhEleadingPt; //! - - TH1F* fhMinRegPtDist; - TH1F* fhRegionMultMin; - TH1F* fhMinRegAvePt; - TH1F* fhMinRegSumPt; - TH1F* fhMinRegMaxPtPart; - TH1F* fhMinRegSumPtvsMult; - - TH1F* fhdNdEta_PhiDist; //! - TH2F* fhFullRegPartPtDistVsEt; //! - TH2F* fhTransRegPartPtDistVsEt; //! - - TH1F* fhRegionSumPtMaxVsEt; //! - TH1I* fhRegionMultMax; //! - TH1F* fhRegionMultMaxVsEt; //! - TH1F* fhRegionSumPtMinVsEt; //! - TH1F* fhRegionMultMinVsEt; //! - TH1F* fhRegionAveSumPtVsEt; //! - TH1F* fhRegionDiffSumPtVsEt; //! - - TH1F* fhRegionAvePartPtMaxVsEt; //! - TH1F* fhRegionAvePartPtMinVsEt; //! - TH1F* fhRegionMaxPartPtMaxVsEt; //! - - // TH2F* fhValidRegion; //! test to be canceled - - ClassDef( AliAnalysisTaskUE, 1) // Analysis task for Underlying Event analysis -}; - #endif + +