Initial structure for PWG4JetTasks (for PWG4 Task Force) incl. steering macro and...
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Jul 2008 11:11:54 +0000 (11:11 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Jul 2008 11:11:54 +0000 (11:11 +0000)
PWG4/AliAnalysisTaskUE.cxx [new file with mode: 0644]
PWG4/AliAnalysisTaskUE.h [new file with mode: 0644]
PWG4/Makefile
PWG4/PROOF-INF.PWG4JetTasks/BUILD.sh [new file with mode: 0755]
PWG4/PROOF-INF.PWG4JetTasks/SETUP.C [new file with mode: 0644]
PWG4/PWG4JetTasksLinkDef.h [new file with mode: 0644]
PWG4/libPWG4JetTasks.pkg [new file with mode: 0644]
PWG4/macros/AnalysisTrainCAF.C [new file with mode: 0644]

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