Add data members for jet background, enable filling (Leticia Mendez)
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 4 Jul 2010 18:33:21 +0000 (18:33 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 4 Jul 2010 18:33:21 +0000 (18:33 +0000)
JETAN/AliAnalysisTaskJets.cxx
JETAN/AliFastJetFinder.cxx
JETAN/AliJetBkg.cxx
JETAN/AliJetBkg.h

index dd48fe8..8bb2fb5 100644 (file)
@@ -153,6 +153,7 @@ void AliAnalysisTaskJets::UserCreateOutputObjects()
        if(fAODExtension)fJetFinder->ConnectAODNonStd(fAODExtension->GetAOD(), fNonStdBranch.Data()); 
       }
       else{
+       if (fDebug > 1) printf("Connecting Non Std Branch AOD %p %s \n",AODEvent(),fNonStdBranch.Data());
        fJetFinder->ConnectAODNonStd(AODEvent(), fNonStdBranch.Data()); 
       }
     }
index 14a1eb9..434b989 100644 (file)
@@ -427,6 +427,18 @@ Bool_t AliFastJetFinder::ProcessEvent()
 
   fJetBkg->SetHeader(fHeader);
   fJetBkg->SetReader(fReader);
+  Double_t sigma1,meanarea1,sigma2,meanarea2;
+  Double_t bkg1,bkg2;
+  fJetBkg->SetFastJetInput(fInputFJ);
+  fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1);
+  fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2);
+  fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1);
+  fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2);
+
+
+
+
+
   /*
   fJetBkg->SetFastJetInput(fInputFJ);
   Double_t bkg1=fJetBkg->BkgFastJet();
@@ -469,6 +481,15 @@ Bool_t AliFastJetFinder::ProcessEvent2()
   fJetBkg->SetHeader(fHeader);
   fJetBkg->SetReader(fReader);
   fJetBkg->SetFastJetInput(fInputFJ);
+  Double_t sigma1,meanarea1,sigma2,meanarea2;
+  Double_t bkg1,bkg2;
+  fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1);
+  fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2);
+  fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1);
+  fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2);
+
+
+
 //  Double_t bkg1=fJetBkg->BkgFastJet();
 //  Double_t bkg2=fJetBkg->BkgChargedFastJet();
 //  Double_t bkg3=fJetBkg->BkgFastJetCone(fAODjets);
index 27d5b89..76638d8 100644 (file)
@@ -83,6 +83,211 @@ AliJetBkg& AliJetBkg::operator=(const AliJetBkg& source){
     return *this;
 
 }
+
+
+//___________________________________________________________________
+  void AliJetBkg::BkgFastJetb(Double_t& rho,Double_t& sigma, 
+Double_t& meanarea){
+  
+  cout<<"===============  AliJetBkg::BkgFastJetb()  =========== "<<endl;
+  vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
+  
+  //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
+  
+
+   
+  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+  double rParamBkg = header->GetRparamBkg(); //Radius for background calculation
+
+  Double_t medianb,sigmab,meanareab;
+  CalcRhob(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All");
+  rho=medianb;
+  sigma=sigmab;
+  meanarea=meanareab;
+}
+//_________________________________________________________________
+
+  void AliJetBkg::BkgFastJetWoHardest(Double_t& rho,Double_t& sigma, 
+Double_t& meanarea){
+  
+  cout<<"===============  AliJetBkg::BkgWoHardest()  =========== "<<endl;
+  vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
+  
+  //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
+  
+
+   
+  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+  double rParamBkg = header->GetRparamBkg(); //Radius for background calculation  
+  Double_t medianb,sigmab,meanareab;
+  CalcRhoWoHardest(medianb,sigmab,meanareab,inputParticles,rParamBkg,"All");
+  rho=medianb;
+  sigma=sigmab;
+  meanarea=meanareab;
+
+  }
+
+//____________________________________________________________________
+
+ void AliJetBkg::CalcRhob(Double_t& median,Double_t& 
+sigma,Double_t& 
+meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t 
+rParamBkg,TString method){
+  //calculate rho using the fastjet method
+
+  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+  Bool_t debug  = header->GetDebug();     // debug option
+
+  fastjet::Strategy strategy = header->GetStrategy();
+  fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
+  fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); 
+  fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
+
+  // create an object that specifies how we to define the area
+  fastjet::AreaDefinition areaDef;
+  double ghostEtamax = header->GetGhostEtaMax(); 
+  double ghostArea   = header->GetGhostArea(); 
+  int    activeAreaRepeats = header->GetActiveAreaRepeats(); 
+
+  // now create the object that holds info about ghosts
+
+  
+
+  fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
+  // and from that get an area definition
+  fastjet::AreaType areaType = header->GetAreaType();
+  areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
+  
+  //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
+  fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef,areaDef);
+  TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
+  comment+= "Jet definition: ";
+  comment+= TString(jetDef.description());
+  //  comment+= ". Area definition: ";
+  //comment+= TString(areaDef.description());
+  comment+= ". Strategy adopted by FastJet: ";
+  comment+= TString(clust_seq.strategy_string());
+  header->SetComment(comment);
+  if(debug){
+    cout << "--------------------------------------------------------" << endl;
+    cout << comment << endl;
+    cout << "--------------------------------------------------------" << endl;
+  }
+
+  vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.);
+  vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); 
+
+  //cout<<"# of BKG jets = "<<jets.size()<<endl;
+      //  for (size_t j = 0; j < jets.size(); j++) { // loop for jets   
+      //printf("BKG Jet found %5d %9.5f %8.5f %10.3f %4.4f 
+      //\n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(),
+      //clust_seq.area(jets[j]));
+      // }
+  
+  
+  double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
+
+    phiMin = 0;
+    phiMax = 2*TMath::Pi();
+  rapMax = ghostEtamax - rParamBkg;
+  rapMin = - ghostEtamax + rParamBkg;
+
+
+  fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
+
+
+  double medianb, sigmab, meanareab;
+  clust_seq.get_median_rho_and_sigma(inclusiveJets, range, false, medianb, sigmab, meanareab, false);
+  median=medianb;
+  sigma=sigmab;
+  meanarea=meanareab; 
+
+ }
+
+//____________________________________________________________________
+
+
+ void AliJetBkg::CalcRhoWoHardest(Double_t& median,Double_t& 
+sigma,Double_t& 
+meanarea,vector<fastjet::PseudoJet> inputParticles,Double_t 
+rParamBkg,TString method){
+  //calculate rho using the fastjet method
+
+  AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
+  Bool_t debug  = header->GetDebug();     // debug option
+
+  fastjet::Strategy strategy = header->GetStrategy();
+  fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
+  fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); 
+  fastjet::JetDefinition jetDef(algorithm, rParamBkg, recombScheme, strategy);
+
+  // create an object that specifies how we to define the area
+  fastjet::AreaDefinition areaDef;
+  double ghostEtamax = header->GetGhostEtaMax(); 
+  double ghostArea   = header->GetGhostArea(); 
+  int    activeAreaRepeats = header->GetActiveAreaRepeats(); 
+
+  // now create the object that holds info about ghosts
+
+
+  fastjet::GhostedAreaSpec ghost_spec(ghostEtamax, activeAreaRepeats, ghostArea);
+  // and from that get an area definition
+  fastjet::AreaType areaType = header->GetAreaType();
+  areaDef = fastjet::AreaDefinition(areaType,ghost_spec);
+  //cout<<"rParamBkg="<<rParamBkg<<"  ghostEtamax="<<ghostEtamax<<"  ghostArea="<<ghostArea<<" areadef="<<TString(areaDef.description())<<endl;
+  //fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef);
+  fastjet::ClusterSequenceArea clust_seq(inputParticles,jetDef,areaDef);
+  TString comment = "Running FastJet algorithm for BKG calculation with the following setup. ";
+  comment+= "Jet definition: ";
+  comment+= TString(jetDef.description());
+  //  comment+= ". Area definition: ";
+  //comment+= TString(areaDef.description());
+  comment+= ". Strategy adopted by FastJet: ";
+  comment+= TString(clust_seq.strategy_string());
+  header->SetComment(comment);
+  if(debug){
+    cout << "--------------------------------------------------------" << endl;
+    cout << comment << endl;
+    cout << "--------------------------------------------------------" << endl;
+  }
+
+  vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(0.);
+  vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); 
+  vector<fastjet::PseudoJet> jets2=sorted_by_pt(inclusiveJets);
+  if(jets2.size()>=2) jets2.erase(jets2.begin(),jets2.begin()+1);
+  
+   
+  
+  double phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
+
+    phiMin = 0;
+    phiMax = 2*TMath::Pi();
+  rapMax = ghostEtamax - rParamBkg;
+  rapMin = - ghostEtamax + rParamBkg;
+
+
+  fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
+
+
+  double medianb, sigmab, meanareab;
+  clust_seq.get_median_rho_and_sigma(jets2, range, false, medianb, sigmab, 
+meanareab, false);
+  median=medianb;
+  sigma=sigmab;
+  meanarea=meanareab; 
+
+  
+ }
+
+
 //___________________________________________________________________
 Float_t AliJetBkg::BkgFastJet(){
   
index 61ff8d0..38a5299 100644 (file)
@@ -31,6 +31,8 @@ class AliJetBkg : public TObject
     void SetHeader(AliJetHeader *header)  {fHeader=header;}
     void SetReader(AliJetReader *reader)  {fReader=reader;}
     void SetFastJetInput(AliFastJetInput *fjinput)  {fInputFJ=fjinput;}
+    void BkgFastJetb(Double_t& x,Double_t& y, Double_t& z);
+    void BkgFastJetWoHardest(Double_t& x,Double_t& y, Double_t& z);
     Float_t BkgFastJet();
     Float_t BkgChargedFastJet();
     Float_t BkgStat();
@@ -43,7 +45,12 @@ class AliJetBkg : public TObject
     
  private:
     Double_t CalcRho(vector<fastjet::PseudoJet> input_particles,Double_t RparamBkg,TString method);
-
+    void CalcRhob(Double_t& median, Double_t& sigma, Double_t& 
+meanarea,vector<fastjet::PseudoJet> input_particles,Double_t RparamBkg,TString 
+method);
+    void CalcRhoWoHardest(Double_t& median, Double_t& sigma, Double_t& 
+meanarea,vector<fastjet::PseudoJet> input_particles,Double_t RparamBkg,TString 
+method);
     AliJetReader *fReader;   //! reader
     AliJetHeader *fHeader;   //! header
     AliFastJetInput *fInputFJ; //! input particles