]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - JETAN/AliJetBkg.cxx
TRD changes for the makeOCDB.C
[u/mrichter/AliRoot.git] / JETAN / AliJetBkg.cxx
index 27d5b894bc0d9d12dacd0133932f7703ba8a92bf..ad3c944aee10746b3f0ca753b3a151cdbd8e601c 100644 (file)
@@ -83,9 +83,215 @@ AliJetBkg& AliJetBkg::operator=(const AliJetBkg& source){
     return *this;
 
 }
+
+
 //___________________________________________________________________
-Float_t AliJetBkg::BkgFastJet(){
+  void AliJetBkg::BkgFastJetb(Double_t& rho,Double_t& sigma, 
+Double_t& meanarea){
+
+    AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; 
+    Bool_t debug  = header->GetDebug();     // debug option 
+    if(debug)cout<<"===============  AliJetBkg::BkgFastJetb()  =========== "<<endl;
+  vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
+  
+  //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;
+  
+
+   
+
+  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){
+    AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; 
+    Bool_t debug  = header->GetDebug();     // debug option 
+    if(debug)cout<<"===============  AliJetBkg::BkgWoHardest()  =========== "<<endl;
+  vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
+  
+  //cout<<"printing inputParticles for BKG "<<inputParticles.size()<<endl;   
+  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());
+  comment+= Form("Method: %s",method.Data());
+  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());
+  comment+= Form("Method: %s",method.Data());
+  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(){
+// Return background  
   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
   Bool_t debug  = header->GetDebug();     // debug option
 
@@ -107,7 +313,7 @@ Float_t AliJetBkg::BkgFastJet(){
 }
 //___________________________________________________________________
 Float_t  AliJetBkg::BkgChargedFastJet(){
-
+// Background for charged jets
   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
   Bool_t debug  = header->GetDebug();     // debug option
 
@@ -126,7 +332,7 @@ Float_t  AliJetBkg::BkgChargedFastJet(){
 
   Double_t rho=CalcRho(inputParticlesCharged,rParam,"Charg");
 
-  cout<<"-------- rho (from CHARGED part)="<<rho<<endl; 
+  if(debug)cout<<"-------- rho (from CHARGED part)="<<rho<<endl; 
   return rho;
 }
 
@@ -191,7 +397,7 @@ Float_t AliJetBkg::BkgFastJetCone(TClonesArray* fAODJets)
   Float_t rCone=0.4;
   
   if(nJ==1) { 
-    AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
+    AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0));
     jeteta=jettmp->Eta();
     jetphi=jettmp->Phi();
     acc=EmcalAcceptance(jeteta,jetphi,rCone);
@@ -201,8 +407,8 @@ Float_t AliJetBkg::BkgFastJetCone(TClonesArray* fAODJets)
 
   
   if(nJ>=2) { 
-    AliAODJet *jettmp = dynamic_cast<AliAODJet*>(fAODJets->At(0));
-    AliAODJet *jettmp1 = dynamic_cast<AliAODJet*>(fAODJets->At(1));
+    AliAODJet *jettmp = (AliAODJet*)(fAODJets->At(0));
+    AliAODJet *jettmp1 = (AliAODJet*)(fAODJets->At(1));
     jeteta=jettmp->Eta();
     jetphi=jettmp->Phi();
     jeteta1=jettmp1->Eta();
@@ -387,8 +593,8 @@ Double_t  AliJetBkg::BkgFunction(Double_t */*x*/,Double_t */*par*/){
 }
 
 
-Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius){
+Bool_t AliJetBkg::EmcalAcceptance(const Float_t eta, const Float_t phi, const Float_t radius) const{
+// Apply emcal acceptance cuts
   Float_t meanPhi=190./180.*TMath::Pi()-110./180.*TMath::Pi()/2;
   Float_t deltaphi=110./180.*TMath::Pi();
   Float_t phicut=deltaphi/2.-radius;