enabling track references for Momentum array (M. Estienne)

diff --git a/JETAN/AliFastJetFinder.cxx b/JETAN/AliFastJetFinder.cxx
index 4d20f9828185c8b57d3191df6afc00222e279389..f15989bffd59ea8ea7adb715a2276608bcf8ac6b 100644 (file)
--- a/JETAN/AliFastJetFinder.cxx
+++ b/JETAN/AliFastJetFinder.cxx
@@ -93,11 +93,10 @@ void AliFastJetFinder::FindJets()
   if(debug)cout<<"----------in AliFastJetFinder::FindJets() ------------------"<<endl;
 
   // check if we are reading AOD jets
-  /* NOT used
   TRefArray *refs = 0;
   Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
   if (fromAod) { refs = fReader->GetReferences(); }
-  */
+
   // RUN ALGORITHM  
   // read input particles -----------------------------
 
@@ -110,6 +109,7 @@ void AliFastJetFinder::FindJets()
   // create an object that represents your choice of jet algorithm, and 
   // the associated parameters
   double rParam = header->GetRparam();
+  double rBkgParam = header->GetRparamBkg();
   fastjet::Strategy strategy = header->GetStrategy();
   fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
   fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); 
@@ -127,41 +127,42 @@ void AliFastJetFinder::FindJets()
   fastjet::AreaType areaType = header->GetAreaType();
   areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
   
+  //***************************** JETS FINDING
+  // run the jet clustering with the above jet definition
+  fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
 
+  vector<fastjet::PseudoJet> jets;
 
-
-
-  if(bgMode) // BG subtraction!!!!!!!!!! Not entering here
+  if(bgMode) // Do BG subtraction directly with the same algorithm (cambridge or kt) for jet signal and background
     {
-      //***************************** JETS FINDING AND EXTRACTION
+      //***************************** CLUSTER JETS FINDING FOR RHO ESTIMATION
       // run the jet clustering with the above jet definition
-      fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
+         fastjet::JetAlgorithm algorithmBkg = header->GetBGAlgorithm();
+         fastjet::JetDefinition jetDefBkg(algorithmBkg, rBkgParam, recombScheme, strategy);
+      fastjet::ClusterSequenceArea clust_seq_bkg(inputParticles, jetDefBkg, areaDef);
 
       // save a comment in the header
-      
       TString comment = "Running FastJet algorithm with the following setup. ";
       comment+= "Jet definition: ";
       comment+= TString(jetDef.description());
+      comment+= "Jet bckg definition: ";
+      comment+= TString(jetDefBkg.description());
       comment+= ". Area definition: ";
       comment+= TString(areaDef.description());
-      comment+= ". Strategy adopted by FastJet: ";
+      comment+= ". Strategy adopted by FastJet and bkg: ";
       comment+= TString(clust_seq.strategy_string());
       header->SetComment(comment);
       if(debug){
-       cout << "--------------------------------------------------------" << endl;
-       cout << comment << endl;
-       cout << "--------------------------------------------------------" << endl;
+          cout << "--------------------------------------------------------" << endl;
+          cout << comment << endl;
+          cout << "--------------------------------------------------------" << endl;
       }
       //header->PrintParameters();
       
-      
       // extract the inclusive jets with pt > ptmin, sorted by pt
       double ptmin = header->GetPtMin(); 
-      vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
-      
-      //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
-      
-      
+      vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets();
+
       //subtract background // ===========================================
       // set the rapididty , phi range within which to study the background 
       double rapMax = header->GetRapMax(); 
@@ -170,121 +171,128 @@ void AliFastJetFinder::FindJets()
       double phiMin = header->GetPhiMin();
       fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
       
-      // subtract background
-      vector<fastjet::PseudoJet> subJets =  clust_seq.subtracted_jets(range,ptmin);  
-      
-      // print out
+      // Extract rho and sigma
+      Double_t rho = 0.;
+      Double_t sigma = 0.;
+      Double_t meanarea = 0.;
+      Bool_t Use4VectorArea = header->Use4VectorArea();
+      vector<fastjet::PseudoJet> bkgJets = clust_seq_bkg.inclusive_jets();
+      clust_seq_bkg.get_median_rho_and_sigma(bkgJets,range, Use4VectorArea, rho, sigma, meanarea, false);
+
+         // subtract background and extract jets bkg subtracted
+      vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(rho,ptmin);
+
+         // print out
       //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n";
       //cout << "---------------------------------------\n";
       //cout << endl;
       //printf(" ijet   rap      phi        Pt         area  +-   err\n");
       
       // sort jets into increasing pt
-      vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);  
-      for (size_t j = 0; j < jets.size(); j++) { // loop for jets
-       
-       double area     = clust_seq.area(jets[j]);
-       double areaError = clust_seq.area_error(jets[j]);
-       
-       if(debug) printf("Jet found %5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n", (Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, areaError);
-       
-       // go to write AOD  info
-       AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
-       //cout << "Printing jet " << endl;
-       if(debug) aodjet.Print("");
-       //cout << "Adding jet ... " ;
-       AddJet(aodjet);
-       //cout << "added \n" << endl;
-       
-      }
+      jets = sorted_by_pt(subJets);  
+
     }
+  else { // No BG subtraction!!!!!!!! Default header is bgmode=0.
 
+         // save a comment in the header
+         TString comment = "Running FastJet algorithm with the following setup. ";
+         comment+= "Jet definition: ";
+         comment+= TString(jetDef.description());
+         comment+= ". Strategy adopted by FastJet: ";
+         comment+= TString(clust_seq.strategy_string());
+         header->SetComment(comment);
+         if(debug){
+                 cout << "--------------------------------------------------------" << endl;
+                 cout << comment << endl;
+                 cout << "--------------------------------------------------------" << endl;
+         }
+         //header->PrintParameters();
 
+         // extract the inclusive jets with pt > ptmin, sorted by pt
+         double ptmin = header->GetPtMin();
+         vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
 
+         jets = sorted_by_pt(inclusiveJets); // Added by me
 
+  }
 
+  for (size_t j = 0; j < jets.size(); j++) { // loop for jets
 
-  else { // No BG subtraction!!!!!!!! Default header is bgmode=0.
-   
-    TClonesArray* fUnit = fReader->GetUnitArray(); //Big mmentum array
-    if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
-    Int_t         nIn = fUnit->GetEntries();
+         double area      = clust_seq.area(jets[j]);
+         double areaError = clust_seq.area_error(jets[j]);
 
- 
-    //fastjet::ClusterSequence clust_seq(inputParticles, jetDef); 
-    fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
-   
-    // save a comment in the header
-    
-    TString comment = "Running FastJet algorithm with the following setup. ";
-    comment+= "Jet definition: ";
-    comment+= TString(jetDef.description());
-    comment+= ". Strategy adopted by FastJet: ";
-    comment+= TString(clust_seq.strategy_string());
-    header->SetComment(comment);
-    if(debug){
-      cout << "--------------------------------------------------------" << endl;
-      cout << comment << endl;
-      cout << "--------------------------------------------------------" << endl;
-    }
-    //header->PrintParameters();
-  
-      // extract the inclusive jets with pt > ptmin, sorted by pt
-    double ptmin = header->GetPtMin(); 
-    vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
-    
-    //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
-    
-    vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); // Added by me
-    for (size_t j = 0; j < jets.size(); j++) { // loop for jets     // Added by me
-      
-      if(debug) printf("Jet found %5d %9.5f %8.5f %10.3f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp());
+         if(debug) printf("Jet found %5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n", (Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, areaError);
 
       vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]);
       int nCon= constituents.size();
       TArrayI ind(nCon);
-      Double_t area=clust_seq.area(jets[j]);
-      // go to write AOD  info
-      AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
-      aodjet.SetEffArea(area,0);
-      //cout << "Printing jet " << endl;
-      if(debug) aodjet.Print("");
-      Int_t count=0;
-      for (int i=0; i < nCon; i++)
-       {
-       fastjet::PseudoJet mPart=constituents[i];
-       ind[i]=mPart.user_index();
-       //      cout<<i<<"  index="<<ind[i]<<endl;
-       
-       //internal loop over all the unit cells
-       Int_t ipart = 0;
-       //      count=0;
-       for(Int_t ii=0; ii<nIn; ii++) 
+      
+      if ((jets[j].eta() > (header->GetJetEtaMax())) ||
+          (jets[j].eta() < (header->GetJetEtaMin())) ||
+          (jets[j].phi() > (header->GetJetPhiMax())) ||
+          (jets[j].phi() < (header->GetJetPhiMin())) ||
+          (jets[j].perp() < header->GetPtMin())) continue; // acceptance eta range and etmin
+
+         // go to write AOD  info
+         AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
+      aodjet.SetEffArea(area,areaError);
+         //cout << "Printing jet " << endl;
+         if(debug) aodjet.Print("");
+
+         Int_t count=0;
+         for (int i=0; i < nCon; i++)
          {
-           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(ii);
-           if(uArray->GetUnitEnergy()>0.){
-             uArray->SetUnitTrackID(ipart);//used to have the index available in AliJEtBkg
-             if(ipart==ind[i]){
-               TRefArray* trackArray = (TRefArray*)uArray->GetUnitTrackRef();
-               Int_t tracksInCell = trackArray->GetEntries();
-               for(int ji = 0; ji < tracksInCell; ji++){
-                 AliAODTrack * track = (AliAODTrack*)trackArray->At(ji);
-                 aodjet.AddTrack(track);
-               }
-
-               count++;
-             }
-             ipart++;
-           }
-         }
-      }
-
-    
-
-      AddJet(aodjet);
-           
-    } // end loop for jets
-  } 
+                 fastjet::PseudoJet mPart=constituents[i];
+                 ind[i]=mPart.user_index();
+                 //    cout<<i<<"  index="<<ind[i]<<endl;
+
+                 if(fromAod)
+                 {
+                         if(fReader->GetReaderHeader()->GetDetector()==0)
+                         {
+                                 for (Int_t iref = 0; iref < refs->GetEntries(); iref++)
+                                 {
+                                         if(iref==ind[i]){
+                                                 AliAODTrack * track = (AliAODTrack*)refs->At(iref);
+                                                 aodjet.AddTrack(track);
+                                         }
+                                 }
+
+                         } else {
+
+                                 TClonesArray* fUnit = fReader->GetUnitArray(); //Big mmentum array
+                                 if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
+                                 Int_t         nIn = fUnit->GetEntries();
+
+                                 //internal loop over all the unit cells
+                                 Int_t ipart = 0;
+
+                                 for(Int_t ii=0; ii<nIn; ii++)
+                                 {
+                                         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(ii);
+                                         if(uArray->GetUnitEnergy()>0.){
+                                                 uArray->SetUnitTrackID(ipart);//used to have the index available in AliJEtBkg
+                                                 if(ipart==ind[i]){
+                                                         TRefArray* trackArray = (TRefArray*)uArray->GetUnitTrackRef();
+                                                         Int_t tracksInCell = trackArray->GetEntries();
+                                                         for(int ji = 0; ji < tracksInCell; ji++){
+                                                                 AliAODTrack * track = (AliAODTrack*)trackArray->At(ji);
+                                                                 aodjet.AddTrack(track);
+                                                         }
+
+                                                         count++;
+                                                 }
+                                                 ipart++;
+                                         }
+                                 }
+                         }
+                 }
+        } 
+
+ AddJet(aodjet);
+
+
+  } // End loop on jets
 
 }
 

diff --git a/JETAN/AliFastJetHeaderV1.cxx b/JETAN/AliFastJetHeaderV1.cxx
index 08c1af5b2b38d6efb0415e3f23c63f71358abada..a7adaba0bde76345d80addb5cbd7563dae82a2fd 100644 (file)
--- a/JETAN/AliFastJetHeaderV1.cxx
+++ b/JETAN/AliFastJetHeaderV1.cxx
@@ -38,6 +38,7 @@ AliFastJetHeaderV1::AliFastJetHeaderV1():
     fRparam(0.4),
     fRparamBkg(0.4),  
     fAlgorithm(fastjet::kt_algorithm),
+    fBGAlgorithm(fastjet::kt_algorithm),
     fStrategy(fastjet::Best),
     fRecombScheme(fastjet::BIpt_scheme),
     fGhostEtaMax(2.0),
@@ -50,7 +51,8 @@ AliFastJetHeaderV1::AliFastJetHeaderV1():
     fRapMin(0.),
     fPhiMax(TMath::TwoPi()),
     fPhiMin(0),
-    fBGMode(0)
+    fBGMode(0),
+    fUse4VectorArea(kFALSE)
 {
   // Constructor
   
@@ -71,6 +73,7 @@ void AliFastJetHeaderV1::PrintParameters() const
   cout << "-- Jet Definition --- " << endl;
   cout << "R " << fRparam << endl;
   cout << "Jet Algorithm " << fAlgorithm << endl; 
+  cout << "Bkg Jet Algorithm " << fBGAlgorithm << endl;
   cout << "Strategy " << fStrategy << endl;  
   cout << "Recombination Scheme " << fRecombScheme << endl; 
   

diff --git a/JETAN/AliFastJetHeaderV1.h b/JETAN/AliFastJetHeaderV1.h
index d8f78596239c9c4041d59e4caee64db2d62ce0dd..8c8ff1e396b44c106c78b14775101442c45265e3 100644 (file)
--- a/JETAN/AliFastJetHeaderV1.h
+++ b/JETAN/AliFastJetHeaderV1.h
@@ -28,6 +28,7 @@ class AliFastJetHeaderV1 : public AliJetHeader
   // Getters
   Double_t                     GetRparam()            const {return fRparam;}
   fastjet::JetAlgorithm        GetAlgorithm()         const {return fAlgorithm;}
+  fastjet::JetAlgorithm        GetBGAlgorithm()       const {return fBGAlgorithm;} 
   fastjet::Strategy            GetStrategy()          const {return fStrategy;}
   fastjet::RecombinationScheme GetRecombScheme()      const {return fRecombScheme;}
   Double_t                     GetGhostEtaMax()       const {return fGhostEtaMax;}
@@ -43,10 +44,12 @@ class AliFastJetHeaderV1 : public AliJetHeader
   Float_t                      GetMinCellEt()         const {return fMinCellEt;} 
   Bool_t                       GetBGMode()            const {return fBGMode;}
   Double_t                     GetRparamBkg()            const {return fRparamBkg;}
+  Bool_t                       Use4VectorArea()       const {return fUse4VectorArea;}  
 
   // Setters
   void SetRparam(Double_t f)                           {fRparam = f;}
   void SetAlgorithm(fastjet::JetAlgorithm f)           {fAlgorithm = f;}
+  void SetBGAlgorithm(fastjet::JetAlgorithm f)         {fBGAlgorithm = f;}
   void SetStrategy(fastjet::Strategy f)                {fStrategy = f;}
   void SetRecombScheme(fastjet::RecombinationScheme f) {fRecombScheme = f;}
   void SetGhostEtaMax(Double_t f)                      {fGhostEtaMax = f;}
@@ -57,7 +60,8 @@ class AliFastJetHeaderV1 : public AliJetHeader
   void SetPhiRange(Double_t fmin, Double_t fmax)       {fPhiMin = fmin; fPhiMax = fmax;}
   void SetPtMin(Double_t ptmin)                        {fPtMin = ptmin;}
   void SetBGMode(Bool_t bgmode)                        {fBGMode = bgmode;}
-  
+  void SetUse4VectorArea()                             {fUse4VectorArea = kTRUE;}
+ 
   void SetComment(TString com) {fComment=com;}
   void SetComment(const char* com) {AliJetHeader::SetComment(com);}
   
@@ -72,6 +76,7 @@ class AliFastJetHeaderV1 : public AliJetHeader
   Double_t fRparam;   // R param
   Double_t fRparamBkg;//R param for bkg calculation
   fastjet::JetAlgorithm fAlgorithm; //fastjet::kt_algorithm
+  fastjet::JetAlgorithm fBGAlgorithm; //fastjet::kt_algorithm
   fastjet::Strategy fStrategy;  //= fastjet::Best;
   fastjet::RecombinationScheme fRecombScheme; // = fastjet::BIpt_scheme;
   
@@ -91,8 +96,9 @@ class AliFastJetHeaderV1 : public AliJetHeader
   Double_t fRapMax, fRapMin; // rapidity range of background sub 
   Double_t fPhiMax, fPhiMin; // phi range of background sub
   Bool_t   fBGMode;          // Do we subtract BG or not?
-  
-  ClassDef(AliFastJetHeaderV1,2)
+  Bool_t   fUse4VectorArea;  // Toggle use of 4-vector area 
+
+  ClassDef(AliFastJetHeaderV1,3)
 };
  
 #endif

diff --git a/JETAN/AliJetHeader.cxx b/JETAN/AliJetHeader.cxx
index c03fb175cc03fed538ec536e52e6e6877e7a626c..22d53f49ac1198188fd210111ec49ded7c82ba15 100755 (executable)
--- a/JETAN/AliJetHeader.cxx
+++ b/JETAN/AliJetHeader.cxx
@@ -18,7 +18,7 @@
 // Stores a comment which describes the jet analysis
 // Author: jgcn@mda.cinvestav.mx
 //---------------------------------------------------------------------
- 
+
 #include "AliJetHeader.h"
 ClassImp(AliJetHeader)
  
@@ -30,6 +30,8 @@ AliJetHeader::AliJetHeader():
   fComment("No comment"),
   fJetEtaMax(0.5),
   fJetEtaMin(-0.5),
+  fJetPhiMax(2*TMath::Pi()),
+  fJetPhiMin(0.),
   fDebug(0)
 {
   //
@@ -45,6 +47,8 @@ AliJetHeader::AliJetHeader(const char * name):
   fComment("No comment"),
   fJetEtaMax(0.5),
   fJetEtaMin(-0.5),
+  fJetPhiMax(2*TMath::Pi()),
+  fJetPhiMin(0.),
   fDebug(0)
 {
   //

diff --git a/JETAN/AliJetHeader.h b/JETAN/AliJetHeader.h
index b3841374c35eb63894f837dea98a03ab2baa3f59..645a643ab973f3033f47c591a79c2b09efcb2469 100755 (executable)
--- a/JETAN/AliJetHeader.h
+++ b/JETAN/AliJetHeader.h
@@ -13,6 +13,7 @@
  
 #include <TNamed.h>
 #include <TString.h>
+#include <TMath.h>
  
 class AliJetHeader : public TNamed
 {
@@ -26,12 +27,16 @@ class AliJetHeader : public TNamed
   virtual TString GetComment() const {return fComment;} 
   virtual Float_t GetJetEtaMax() const {return fJetEtaMax;}
   virtual Float_t GetJetEtaMin() const {return fJetEtaMin;}
+  virtual Float_t GetJetPhiMax() const {return fJetPhiMax;}
+  virtual Float_t GetJetPhiMin() const {return fJetPhiMin;}
   virtual Int_t   GetDebug() const     {return fDebug;    }
 
   // Setters
   virtual void SetComment(const char* com) {fComment=TString(com);}
   virtual void SetJetEtaMax(Float_t eta= 0.5) {fJetEtaMax = eta;}
   virtual void SetJetEtaMin(Float_t eta= -0.5) {fJetEtaMin = eta;}
+  virtual void SetJetPhiMax(Float_t phi= 2*TMath::Pi()) {fJetPhiMax = phi;}
+  virtual void SetJetPhiMin(Float_t phi= 0.) {fJetPhiMin = phi;}
   virtual void SetDebug(Int_t debug) {fDebug = debug;}
   
 
@@ -41,9 +46,11 @@ protected:
   TString fComment;         // a comment 
   Float_t fJetEtaMax;       // maximum eta for the jet
   Float_t fJetEtaMin;       // minimum eta for the jet
+  Float_t fJetPhiMax;       // maximum phi for the jet
+  Float_t fJetPhiMin;       // minimum phi for the jet
   Int_t   fDebug;           // debug option
   
-  ClassDef(AliJetHeader,2)
+  ClassDef(AliJetHeader,3)
 };
  
 #endif