update GetJetMasS() function

clean-up

clean-up

diff --git a/PWGJE/EMCALJetTasks/AliJetShape.h b/PWGJE/EMCALJetTasks/AliJetShape.h
index 8fe1e79a311452eb863d924f5fc48c8bc5a000e8..c263db9eb49ae5623b77ec14556351ec897e8913 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliJetShape.h
+++ b/PWGJE/EMCALJetTasks/AliJetShape.h
@@ -14,16 +14,6 @@
 
 using namespace std;
 
-/*
-const Int_t knbins = 130;
-const Double_t kxmin = -0.8;
-const Double_t kxmax = 0.5;
-static Int_t GetBin(Double_t x) {Double_t mdx = (kxmax-kxmin)/(double)knbins; Int_t bin = TMath::FloorNint((x-kxmin)/mdx); return bin;}
-static Double_t gfcn[knbins] = {0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000001, 0.000005, 0.000020, 0.000069, 0.000226, 0.000696, 0.002014, 0.005473, 0.013977, 0.033528, 0.075556, 0.159953, 0.318105, 0.594298, 1.043025, 1.719657, 2.663457, 3.875307, 5.296916, 6.801375, 8.204024, 9.296377, 9.895942, 9.895942, 9.296377, 8.204024, 6.801375, 5.296916, 3.875307, 2.663457, 1.719657, 1.043025, 0.594298, 0.318105, 0.159953, 0.075556, 0.033528, 0.013977, 0.005473, 0.002014, 0.000696, 0.000226, 0.000069, 0.000020, 0.000005, 0.000001, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000};
-
-static Double_t efcn[knbins] = {0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000001, 0.000002, 0.000006, 0.000019, 0.000053, 0.000144, 0.000369, 0.000889, 0.002020, 0.004332, 0.008774, 0.016793, 0.030396, 0.052081, 0.084566, 0.130295, 0.190787, 0.265986, 0.353830, 0.450262, 0.549738, 0.646170, 0.734014, 0.809213, 0.869705, 0.915434, 0.947919, 0.969604, 0.983207, 0.991226, 0.995668, 0.997980, 0.999111, 0.999631, 0.999856, 0.999947, 0.999981, 0.999994, 0.999998, 0.999999, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000};
-*/
-
 #ifdef FASTJET_VERSION
 class AliJetShapeMass : public fastjet::FunctionOfPseudoJet<Double32_t>
 {
@@ -40,10 +30,9 @@ class AliJetShapeGRNum : public fastjet::FunctionOfPseudoJet<Double32_t>
   virtual std::string description() const{return "Numerator angular structure function";}
   //  static Int_t GetBin(Double_t x) {Int_t bin = TMath::FloorNint((x-kxmin)/mdx); return bin;}
   Double32_t result(const fastjet::PseudoJet &jet) const {
-    if (!jet.has_constituents()) {
-      Printf("Angular structure can only be applied on jets for which the constituents are known.");
+    if (!jet.has_constituents())
       return 0; //AliFatal("Angular structure can only be applied on jets for which the constituents are known.");
-    }
+
   Double_t A = 0.;
   std::vector<fastjet::PseudoJet> constits = jet.constituents();
   for(UInt_t ic = 0; ic < constits.size(); ++ic) {
@@ -61,15 +50,8 @@ class AliJetShapeGRNum : public fastjet::FunctionOfPseudoJet<Double32_t>
       if(dr2>0.) {
        Double_t dr = TMath::Sqrt(dr2);
        Double_t x = fR-dr;
-       //      Printf("Calculating numerator between %d-%d for r: %f dr: %f",ic,jc,fR,dr);
        //noisy function
        Double_t noise = TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep);
-       /* Int_t bin = GetBin(x); */
-       /* Double_t noise = 0; */
-       /* if(bin<0) noise = gfcn[0]; */
-       /* else if(bin>=knbins) noise = gfcn[knbins-1]; */
-       /* else  noise = gfcn[bin]; */
-       /* Printf("fR: %f dr: %f x=%f bin=%d noise: %f  %f",fR,dr,x,bin,noise,TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep)); */
        A += constits[ic].perp()*constits[jc].perp()*dr2*noise;
       }
     }
@@ -111,11 +93,6 @@ class AliJetShapeGRDen : public fastjet::FunctionOfPseudoJet<Double32_t>
        Double_t x = fR-dr;
        //error function
        Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
-       /* Int_t bin = GetBin(x); */
-       /* Double_t erf = 0; */
-       /* if(bin<0) erf = efcn[0]; */
-       /* else if(bin>=knbins) erf = efcn[knbins-1]; */
-       /* else  erf = efcn[bin]; */
        A += constits[ic].perp()*constits[jc].perp()*dr2*erf;
       }
     }

diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMass.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMass.cxx
index 6bbeec01f87e54209568174b724a40c8eb27fc57..eeec5c89cb69dd2e37ad9a74c9a0765304c28cc4 100644 (file)
--- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMass.cxx
+++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMass.cxx
@@ -40,7 +40,7 @@ AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass() :
   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMass", kTRUE),
   fContainerBase(0),
   fMinFractionShared(0),
-  fJetMassAvg(0),
+  fJetMassType(kRaw),
   fh2PtJet1VsLeadPtAllSel(0),
   fh2PtJet1VsLeadPtTagged(0),
   fh2PtJet1VsLeadPtTaggedMatch(0),
@@ -104,7 +104,7 @@ AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass(const char *name) :
   AliAnalysisTaskEmcalJet(name, kTRUE),  
   fContainerBase(0),
   fMinFractionShared(0),
-  fJetMassAvg(0),
+  fJetMassType(kRaw),
   fh2PtJet1VsLeadPtAllSel(0),
   fh2PtJet1VsLeadPtTagged(0),
   fh2PtJet1VsLeadPtTaggedMatch(0),
@@ -124,7 +124,7 @@ AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass(const char *name) :
 {
   // Standard constructor.
 
- fh2PtJet1VsLeadPtAllSel        = new TH2F*[fNcentBins];
+  fh2PtJet1VsLeadPtAllSel        = new TH2F*[fNcentBins];
   fh2PtJet1VsLeadPtTagged        = new TH2F*[fNcentBins];
   fh2PtJet1VsLeadPtTaggedMatch   = new TH2F*[fNcentBins];
   fh2PtVsMassJet1All             = new TH2F*[fNcentBins];
@@ -355,32 +355,12 @@ Bool_t AliAnalysisTaskEmcalJetMass::FillHistograms()
 //________________________________________________________________________
 Double_t AliAnalysisTaskEmcalJetMass::GetJetMass(AliEmcalJet *jet) {
   //calc subtracted jet mass
-  TLorentzVector vpS = GetSubtractedVector(jet);
-  return vpS.M();
-}
-
-
-//________________________________________________________________________
-TLorentzVector AliAnalysisTaskEmcalJetMass::GetSubtractedVector(AliEmcalJet *jet) {
-  //get subtracted vector
-  TLorentzVector vpS;
-  AliJetContainer *jetCont = GetJetContainer(fContainerBase);
-  TLorentzVector vpBkg = GetBkgVector(jet,jetCont);
-  vpS.SetPxPyPzE(jet->Px()-vpBkg.Px(),jet->Py()-vpBkg.Py(),jet->Pz()-vpBkg.Pz(),jet->E()-vpBkg.E());
-  return vpS;
-}
-
-//________________________________________________________________________
-TLorentzVector AliAnalysisTaskEmcalJetMass::GetBkgVector(AliEmcalJet *jet, AliJetContainer *cont) {
-  //get background vector
-
-  Double_t rho  = cont->GetRhoVal();
-  Double_t rhom = cont->GetRhoMassVal();
-  TLorentzVector vpB;
-  Double_t aphi = jet->AreaPhi();
-  Double_t aeta = jet->AreaEta();
-  vpB.SetPxPyPzE(rho*TMath::Cos(aphi)*jet->Area(),rho*TMath::Sin(aphi)*jet->Area(),(rho+rhom)*TMath::SinH(aeta)*jet->Area(),(rho+rhom)*TMath::CosH(aeta)*jet->Area());
-  return vpB;
+  if(fJetMassType==kRaw)
+    return jet->M();
+  else if(fJetMassType==kDeriv)
+    return jet->GetSecondOrderSubtracted();
+  
+  return 0;
 }
 
 //________________________________________________________________________

diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMass.h b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMass.h
index 09063e2ed5c969d8b4231a269b1018360c541a2d..ec481540490b99571887fd060c37298c3b061a7e 100644 (file)
--- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMass.h
+++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMass.h
@@ -15,6 +15,10 @@ class AliJetContainer;
 
 class AliAnalysisTaskEmcalJetMass : public AliAnalysisTaskEmcalJet {
  public:
+  enum JetMassType {
+    kRaw   = 0,  //mass form anti-kt 4-vector
+    kDeriv = 1   //area based subtracted jet mass
+  };
 
   AliAnalysisTaskEmcalJetMass();
   AliAnalysisTaskEmcalJetMass(const char *name);
@@ -26,7 +30,7 @@ class AliAnalysisTaskEmcalJetMass : public AliAnalysisTaskEmcalJet {
   //Setters
   void SetJetContainerBase(Int_t c)                             { fContainerBase     = c   ; }
   void SetMinFractionShared(Double_t f)                         { fMinFractionShared = f   ; }
-  void SetJetMassAverage(Double_t m)                            { fJetMassAvg        = m   ; }
+  void SetJetMassType(JetMassType t)                            { fJetMassType       = t   ; }
 
  protected:
   Bool_t                              RetrieveEventObjects();
@@ -34,12 +38,10 @@ class AliAnalysisTaskEmcalJetMass : public AliAnalysisTaskEmcalJet {
   Bool_t                              FillHistograms();
 
   Double_t                            GetJetMass(AliEmcalJet *jet);
-  TLorentzVector                      GetSubtractedVector(AliEmcalJet *jet);
-  TLorentzVector                      GetBkgVector(AliEmcalJet *jet, AliJetContainer *cont);
  
   Int_t                               fContainerBase;              // jets to be analyzed
   Double_t                            fMinFractionShared;          // only fill histos for jets if shared fraction larger than X
-  Double_t                            fJetMassAvg;                 // average jet mass
+  JetMassType                         fJetMassType;                // jet mass type to be used
 
   TH2F            **fh2PtJet1VsLeadPtAllSel;         //!all jets after std selection vs leading track pt
   TH2F            **fh2PtJet1VsLeadPtTagged;         //!tagged jets vs leading track pt
@@ -62,7 +64,7 @@ class AliAnalysisTaskEmcalJetMass : public AliAnalysisTaskEmcalJet {
   AliAnalysisTaskEmcalJetMass(const AliAnalysisTaskEmcalJetMass&);            // not implemented
   AliAnalysisTaskEmcalJetMass &operator=(const AliAnalysisTaskEmcalJetMass&); // not implemented
 
-  ClassDef(AliAnalysisTaskEmcalJetMass, 4)
+  ClassDef(AliAnalysisTaskEmcalJetMass, 5)
 };
 #endif
 

diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeGR.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeGR.cxx
index 9d8a641673ffe546fc8f2742844cf78d5e22e59a..eb82017ee2d0016c479564f779b01465de575080 100644 (file)
--- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeGR.cxx
+++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeGR.cxx
@@ -535,8 +535,6 @@ Bool_t AliAnalysisTaskJetShapeGR::FillHistograms()
        //now get second derivative vs R and do final calculation
        TArrayF num = jet1->GetGRNumeratorSub();
        TArrayF den = jet1->GetGRDenominatorSub();
-       Printf("Got numerator size: nr: %d",nr);
-       cout << "size: " << num.GetSize() << endl;
        if(num.GetSize()>0) {
          for(Int_t i = 0; i<nr; i++) {
            Double_t r = i*wr + 0.5*wr;
@@ -577,10 +575,8 @@ Bool_t AliAnalysisTaskJetShapeGR::FillTrueJets() {
 
   AliEmcalJet* jet1 = NULL;
   AliJetContainer *jetCont = GetJetContainer(fContainerTrue);
-  if(!jetCont) {
-    Printf("cannot find container");
+  if(!jetCont)
     return kFALSE;
-  }
 
   AliDebug(11,Form("NJets True: %d",jetCont->GetNJets()));
 
@@ -625,7 +621,6 @@ Double_t AliAnalysisTaskJetShapeGR::CalcDeltaGR(AliEmcalJet *jet, Int_t ic, TArr
   AliVParticle *vp2 = 0x0;
   Double_t A = 0.; Double_t B = 0.;
   if(jet->GetNumberOfTracks()<2) return 0.;
-  //  Printf("jet pt: %f  nconst: %d",jet->Pt(),jet->GetNumberOfTracks());
   for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
     vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
     if(!vp1) continue;
@@ -634,7 +629,6 @@ Double_t AliAnalysisTaskJetShapeGR::CalcDeltaGR(AliEmcalJet *jet, Int_t ic, TArr
       if(!vp2) continue;
       Double_t dphi = GetDeltaPhi(vp1->Phi(),vp2->Phi());
       Double_t dr2 = (vp1->Eta()-vp2->Eta())*(vp1->Eta()-vp2->Eta()) + dphi*dphi;
-      //      Printf("dr2: %f  dEta: %f  dphi: %f",dr2,vp1->Eta()-vp2->Eta(),dphi);
       if(dr2>0.) {
        for(Int_t k = 0; k<nr; k++) {
          Double_t r = k*fDRStep + 0.5*fDRStep;
@@ -670,11 +664,9 @@ Double_t AliAnalysisTaskJetShapeGR::CalcGR(AliEmcalJet *jet, Int_t ic) {
   Double_t wr = 0.04;
   const Int_t nr = TMath::CeilNint(jetCont->GetJetRadius()/wr);
   Double_t grArr[nr];
-  for(Int_t i = 0; i<nr; i++) {
+  for(Int_t i = 0; i<nr; i++)
     grArr[i] = 0.;
-    //Printf("bin up edge %d=%f",i,wr+i*wr);
-  }
-  //  Printf("jet pt: %f  nconst: %d",jet->Pt(),jet->GetNumberOfTracks());
+
   for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
     vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
     for(Int_t j=i; j<jet->GetNumberOfTracks(); j++) {