Adding the background calculation to the cluster task (Leticia)
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 28 Sep 2010 07:35:25 +0000 (07:35 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 28 Sep 2010 07:35:25 +0000 (07:35 +0000)
JETAN/AliAnalysisTaskJetCluster.cxx
JETAN/AliAnalysisTaskJetCluster.h
PWG4/macros/AddTaskJetCluster.C

index 36cf719..161eced 100644 (file)
@@ -55,7 +55,7 @@
 #include "AliJetKineReaderHeader.h"
 #include "AliGenCocktailEventHeader.h"
 #include "AliInputEventHandler.h"
-
+#include "AliAODJetEventBackground.h"
 
 #include "fastjet/PseudoJet.hh"
 #include "fastjet/ClusterSequenceArea.hh"
@@ -94,6 +94,9 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
   fStrategy(fastjet::Best),
   fRecombScheme(fastjet::BIpt_scheme),
   fAreaType(fastjet::active_area), 
+  fGhostArea(0.01),
+  fActiveAreaRepeats(1),
+  fGhostEtamax(1.5),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
@@ -171,6 +174,9 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
   fStrategy(fastjet::Best),
   fRecombScheme(fastjet::BIpt_scheme),
   fAreaType(fastjet::active_area), 
+  fGhostArea(0.01),
+  fActiveAreaRepeats(1),
+  fGhostEtamax(1.5),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
@@ -250,7 +256,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
 
   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 
-
+  
 
   if(fNonStdBranch.Length()!=0)
     {
@@ -262,6 +268,20 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
       TClonesArray *tca = new TClonesArray("AliAODJet", 0);
       tca->SetName(fNonStdBranch.Data());
       AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
+
+   
+      TClonesArray *tcaran = new TClonesArray("AliAODJet", 0);
+      tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
+      AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data());
+
+
+     if(!AODEvent() || !AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
+       AliAODJetEventBackground* evBkg = new AliAODJetEventBackground();
+       evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
+       AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data());
+       
+     }
+
       if(fNonStdFile.Length()!=0){
        // 
        // case that we have an AOD extension we need to fetch the jets from the extended output
@@ -509,12 +529,24 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
   // handle and reset the output jet branch 
   // only need this once
   TClonesArray* jarray = 0;  
+  TClonesArray* jarrayran = 0;  
+  AliAODJetEventBackground*  evBkg = 0;
   if(fNonStdBranch.Length()!=0) {
     if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
     if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()));
     if(jarray)jarray->Delete();    // this is our responsibility, clear before filling again
+    if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
+    if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
+    if(jarrayran)jarrayran->Delete();    // this is our responsibility, clear before filling again
+
+    if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
+    if(!evBkg)  evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
+    if(evBkg)evBkg->Reset(); 
+
   }
 
+
   //
   // Execute analysis for current event
   //
@@ -630,19 +662,25 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
   // run fast jet
   // employ setters for these...
 
-  double ghostEtamax = 0.9;
-  double ghostArea   = 0.01;
-  int    activeAreaRepeats = 1;
   // now create the object that holds info about ghosts                         
-  fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea)\
-;
+  fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
   fastjet::AreaType areaType =   fastjet::active_area;
   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
-  
-  fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
+    fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
   vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
   
+  //range where to compute background
+  Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
+  phiMin = 0;
+  phiMax = 2*TMath::Pi();
+  rapMax = fGhostEtamax - fRparam;
+  rapMin = - fGhostEtamax + fRparam;
+  fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
+
+
   inclusiveJets = clustSeq.inclusive_jets();
   sortedJets    = sorted_by_pt(inclusiveJets);
  
@@ -759,6 +797,27 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
    }
    delete recIter;
+
+      //background estimates:all bckg jets(0) & wo the 2 hardest(1)
+   if(evBkg){
+     vector<fastjet::PseudoJet> jets2=sortedJets;
+     if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); 
+     Double_t bkg1=0;
+     Double_t sigma1=0.;
+     Double_t meanarea1=0.;
+     Double_t bkg2=0;
+     Double_t sigma2=0.;
+     Double_t meanarea2=0.;
+         
+     clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, false);
+     evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
+     clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, false);
+     evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
+   }
+     
+
+
   }
  
  // fill track information
@@ -815,7 +874,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 
  // find the random jets
  vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
- fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef);
+ fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
   
  inclusiveJetsRan = clustSeqRan.inclusive_jets();
  sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
@@ -864,8 +923,8 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
       }  
 
-
-
+    Int_t nAodOutJetsRan = 0;
+     AliAODJet *aodOutJetRan = 0;
     for(int j = 0; j < nRecRan;j++){
       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
       Float_t tmpPt = tmpRec.Pt();
@@ -876,6 +935,27 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
       fh2NConstPtRan->Fill(tmpPt,constituents.size());
       fh2PtNchRan->Fill(nCh,tmpPt);
       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
+
+
+     if(tmpPt>fJetOutputMinPt&&jarrayran){
+       aodOutJetRan =  new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
+       Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
+       
+       aodOutJetRan->SetEffArea(arearan,0);    }
+
+
+
+     
+     for(unsigned int ic = 0; ic < constituents.size();ic++){
+       // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
+       // fh1PtJetConstRec->Fill(part->Pt());
+       if(aodOutJetRan){
+        aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
+       }}
+      
+
+
+
       // correlation
       Float_t tmpPhi =  tmpRec.Phi();
       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
@@ -887,6 +967,18 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
        continue;
       }
     }  
+
+     
+    if(evBkg){
+     Double_t bkg3=0.;
+     Double_t sigma3=0.;
+     Double_t meanarea3=0.;
+     clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, false);
+     evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
+    }
+
+
+
  }
  
 
index 70cb809..4c8722c 100644 (file)
@@ -23,7 +23,8 @@ class AliAODExtension;
 class AliAODJet;
 class AliGenPythiaEventHeader;
 class AliCFManager;
-
+class AliAODJetEventBackground;
+class AliJetFinder;
 class TList;
 class TChain;
 class TH2F;
@@ -47,6 +48,8 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
     virtual void Terminate(Option_t *option);
     virtual Bool_t Notify();
 
+    
+
     virtual void SetUseGlobalSelection(Bool_t b){fUseGlobalSelection = b;}
     virtual void SetAODTrackInput(Bool_t b){fUseAODTrackInput = b;}
     virtual void SetAODMCInput(Bool_t b){fUseAODMCInput = b;}
@@ -55,7 +58,7 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
     virtual void SetTrackTypeRec(Int_t i){fTrackTypeRec = i;}
     virtual void SetTrackPtCut(Float_t x){fTrackPtCut = x;}
     virtual void SetFilterMask(UInt_t i){fFilterMask = i;}
-
+    
     virtual void SetNSkipLeadingRan(Int_t x){fNSkipLeadingRan = x;}
 
     virtual void SetJetOutputBranch(const char *c){fNonStdBranch = c;}
@@ -73,6 +76,11 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
     void SetStrategy(fastjet::Strategy f)                {fStrategy = f;}
     void SetRecombScheme(fastjet::RecombinationScheme f) {fRecombScheme = f;}
     void SetAreaType(fastjet::AreaType f)                {fAreaType = f;}
+    void SetGhostArea(Double_t f) {fGhostArea = f;}
+    void SetActiveAreaRepeats(Int_t f) {fActiveAreaRepeats = f;}
+    void SetGhostEtamax(Double_t f) {fGhostEtamax = f;}
+
+
 
     // Helper
     //
@@ -121,64 +129,66 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
     fastjet::Strategy fStrategy;  //= fastjet::Best;
     fastjet::RecombinationScheme fRecombScheme; // = fastjet::BIpt_scheme;
     fastjet::AreaType fAreaType; 
-
-    TProfile*     fh1Xsec;   // pythia cross section and trials
-    TH1F*         fh1Trials; // trials are added
-    TH1F*         fh1PtHard;  // Pt har of the event...       
-    TH1F*         fh1PtHardNoW;  // Pt har of the event without weigt       
-    TH1F*         fh1PtHardTrials;  // Number of trials 
-
-    TH1F*         fh1NJetsRec; // number of reconstructed jets
-    TH1F*         fh1NConstRec;// number of constiutens in leading jet
-    TH1F*         fh1NConstLeadingRec;// number of constiutens in leading jet
-    TH1F*         fh1PtJetsRecIn;  // Jet pt for all jets
-    TH1F*         fh1PtJetsLeadingRecIn;  // Jet pt for all jets
-    TH1F*         fh1PtJetConstRec;// pt of constituents
+    Double_t fGhostArea;
+    Int_t fActiveAreaRepeats;
+    Double_t fGhostEtamax;
+    TProfile*     fh1Xsec;   //! pythia cross section and trials
+    TH1F*         fh1Trials; //! trials are added
+    TH1F*         fh1PtHard;  //! Pt har of the event...       
+    TH1F*         fh1PtHardNoW;  //! Pt har of the event without weigt       
+    TH1F*         fh1PtHardTrials;  //! Number of trials 
+
+    TH1F*         fh1NJetsRec; //! number of reconstructed jets
+    TH1F*         fh1NConstRec;//! number of constiutens in leading jet
+    TH1F*         fh1NConstLeadingRec;//! number of constiutens in leading jet
+    TH1F*         fh1PtJetsRecIn;  //! Jet pt for all jets
+    TH1F*         fh1PtJetsLeadingRecIn;  //! Jet pt for all jets
+    TH1F*         fh1PtJetConstRec;//! pt of constituents
     TH1F*         fh1PtJetConstLeadingRec;// pt of constituents
-    TH1F*         fh1PtTracksRecIn;  // track pt for all tracks
-    TH1F*         fh1PtTracksLeadingRecIn;  // track pt for all tracks
+    TH1F*         fh1PtTracksRecIn;  //! track pt for all tracks
+    TH1F*         fh1PtTracksLeadingRecIn;  //! track pt for all tracks
 
     // Randomized track histos
-    TH1F*         fh1NJetsRecRan; // number of reconstructed jets from randomized
-    TH1F*         fh1NConstRecRan;// number of constiutens in leading jet
-    TH1F*         fh1PtJetsLeadingRecInRan;  // Jet pt for all jets
-    TH1F*         fh1NConstLeadingRecRan;// number of constiutens in leading jet
-    TH1F*         fh1PtJetsRecInRan;  // Jet pt for all jets
-
-    TH1F*         fh1PtTracksGenIn;  // track pt for all tracks
-    TH1F*         fh1Nch;            // charged particle mult
-
-    TH2F*         fh2NRecJetsPt;            // Number of found jets above threshold
-    TH2F*         fh2NRecTracksPt;          // Number of found tracks above threshold
-    TH2F*         fh2NConstPt;           // number of constituents vs. pt
-    TH2F*         fh2NConstLeadingPt;           // number of constituents vs. pt
-    TH2F*         fh2JetPhiEta;             // jet phi eta
-    TH2F*         fh2LeadingJetPhiEta;      // leading jet phi eta
-    TH2F*         fh2JetEtaPt;              // leading jet eta
-    TH2F*         fh2LeadingJetEtaPt;              // leading jet eta
-    TH2F*         fh2TrackEtaPt;              // track eta
-    TH2F*         fh2LeadingTrackEtaPt;       // leading track eta
-    TH2F*         fh2JetsLeadingPhiEta;     // jet phi eta
-    TH2F*         fh2JetsLeadingPhiPt;      // jet correlation with leading jet
-    TH2F*         fh2TracksLeadingPhiEta;   // track correlation with leading track
-    TH2F*         fh2TracksLeadingPhiPt;    // track correlation with leading track
-    TH2F*         fh2TracksLeadingJetPhiPt; // track correlation with leading Jet
-    TH2F*         fh2JetsLeadingPhiPtW;      // jet correlation with leading jet
-    TH2F*         fh2TracksLeadingPhiPtW;   // track correlation with leading track
-    TH2F*         fh2TracksLeadingJetPhiPtW; // track correlation with leading Jet
-    TH2F*         fh2NRecJetsPtRan;            // Number of found jets above threshold
-    TH2F*         fh2NConstPtRan;           // number of constituents vs. pt
-    TH2F*         fh2NConstLeadingPtRan;           // number of constituents vs. pt
-    TH2F*         fh2PtNch;               // p_T of cluster vs. multiplicity,
-    TH2F*         fh2PtNchRan;            // p_T of cluster vs. multiplicity,random
-    TH2F*         fh2PtNchN;               // p_T of cluster vs. multiplicity, weigthed with constituents
-    TH2F*         fh2PtNchNRan;            // p_T of cluster vs. multiplicity, weigthed with constituents random
-    TH2F*         fh2TracksLeadingJetPhiPtRan; // track correlation with leading Jet
-    TH2F*         fh2TracksLeadingJetPhiPtWRan; // track correlation with leading Jet
+    TH1F*         fh1NJetsRecRan; //! number of reconstructed jets from randomized
+    TH1F*         fh1NConstRecRan;//! number of constiutens in leading jet
+    TH1F*         fh1PtJetsLeadingRecInRan;  //! Jet pt for all jets
+    TH1F*         fh1NConstLeadingRecRan;//! number of constiutens in leading jet
+    TH1F*         fh1PtJetsRecInRan;  //! Jet pt for all jets
+
+    TH1F*         fh1PtTracksGenIn;  //! track pt for all tracks
+    TH1F*         fh1Nch;            //! charged particle mult
+
+    TH2F*         fh2NRecJetsPt;            //! Number of found jets above threshold
+    TH2F*         fh2NRecTracksPt;          //! Number of found tracks above threshold
+    TH2F*         fh2NConstPt;           //! number of constituents vs. pt
+    TH2F*         fh2NConstLeadingPt;           //! number of constituents vs. pt
+    TH2F*         fh2JetPhiEta;             //! jet phi eta
+    TH2F*         fh2LeadingJetPhiEta;      //! leading jet phi eta
+    TH2F*         fh2JetEtaPt;              //! leading jet eta
+    TH2F*         fh2LeadingJetEtaPt;              //! leading jet eta
+    TH2F*         fh2TrackEtaPt;              //! track eta
+    TH2F*         fh2LeadingTrackEtaPt;       //! leading track eta
+    TH2F*         fh2JetsLeadingPhiEta;     //! jet phi eta
+    TH2F*         fh2JetsLeadingPhiPt;      //! jet correlation with leading jet
+    TH2F*         fh2TracksLeadingPhiEta;   //! track correlation with leading track
+    TH2F*         fh2TracksLeadingPhiPt;    //! track correlation with leading track
+    TH2F*         fh2TracksLeadingJetPhiPt; //! track correlation with leading Jet
+    TH2F*         fh2JetsLeadingPhiPtW;      //! jet correlation with leading jet
+    TH2F*         fh2TracksLeadingPhiPtW;   //! track correlation with leading track
+    TH2F*         fh2TracksLeadingJetPhiPtW; //! track correlation with leading Jet
+    TH2F*         fh2NRecJetsPtRan;            //! Number of found jets above threshold
+    TH2F*         fh2NConstPtRan;           //! number of constituents vs. pt
+    TH2F*         fh2NConstLeadingPtRan;           //! number of constituents vs. pt
+    TH2F*         fh2PtNch;               //! p_T of cluster vs. multiplicity,
+    TH2F*         fh2PtNchRan;            //! p_T of cluster vs. multiplicity,random
+    TH2F*         fh2PtNchN;               //! p_T of cluster vs. multiplicity, weigthed with constituents
+    TH2F*         fh2PtNchNRan;            //! p_T of cluster vs. multiplicity, weigthed with constituents random
+    TH2F*         fh2TracksLeadingJetPhiPtRan; //! track correlation with leading Jet
+    TH2F*         fh2TracksLeadingJetPhiPtWRan; //! track correlation with leading Jet
     TList *fHistList; // Output list
    
 
-    ClassDef(AliAnalysisTaskJetCluster, 4) 
+    ClassDef(AliAnalysisTaskJetCluster, 5) 
 };
  
 #endif
index 2638bbf..f65beb5 100644 (file)
@@ -123,7 +123,7 @@ AliAnalysisTaskJetCluster *AddTaskJetCluster(char* bRec,char* bGen ,UInt_t filte
    \r
    if(TMath::Abs(radius-0.4)<0.01){\r
      pwg4spec->SetJetOutputBranch(Form("clusters%s_%s%s",bRec,jf,cRadius));\r
-     pwg4spec->SetJetOutputMinPt(1); // store only jets / clusters above a certain threshold\r
+     pwg4spec->SetJetOutputMinPt(0); // store only jets / clusters above a certain threshold\r
    }\r
 \r
    pwg4spec->SetNSkipLeadingRan(nSkip);\r