first version of B->J/psi->ee analysis (Giuseppe)
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 20 Apr 2007 07:49:45 +0000 (07:49 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 20 Apr 2007 07:49:45 +0000 (07:49 +0000)
PWG3/AliBtoJPSItoEle.cxx [new file with mode: 0644]
PWG3/AliBtoJPSItoEle.h [new file with mode: 0644]
PWG3/AliBtoJPSItoEleAnalysis.cxx [new file with mode: 0644]
PWG3/AliBtoJPSItoEleAnalysis.h [new file with mode: 0644]
PWG3/AliBtoJPSItoEleReco_all.C [new file with mode: 0644]
PWG3/AliPrimJPSItoEleReco_all.C [new file with mode: 0644]

diff --git a/PWG3/AliBtoJPSItoEle.cxx b/PWG3/AliBtoJPSItoEle.cxx
new file mode 100644 (file)
index 0000000..65e5efe
--- /dev/null
@@ -0,0 +1,436 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//----------------------------------------------------------------------------
+//               Implementation of the BtoJPSItoEle class
+//                  for pp and PbPb interactions
+// Note: the two decay tracks are labelled: 0 (positive electron)
+//                                          1 (negative electron)
+//            Origin: G.E. Bruno    giuseppe.bruno@ba.infn.it            
+//  based on Class for charm golden channel (D0->Kpi)
+//----------------------------------------------------------------------------
+
+// #include <Riostream.h> // for debugging
+
+#include <TH1.h>
+#include <TH2.h>
+#include <TCanvas.h>
+#include <TPaveLabel.h>
+#include <TVector3.h>
+
+#include "AliBtoJPSItoEle.h"
+
+ClassImp(AliBtoJPSItoEle)
+
+//----------------------------------------------------------------------------
+AliBtoJPSItoEle::AliBtoJPSItoEle() {
+  // Default constructor
+  
+  fSignal = kFALSE;
+  fJpsiPrimary = kFALSE;
+
+  fEvent = 0;
+
+  fTrkNum[0] = 0;
+  fTrkNum[1] = 0;
+
+  fV1x = 0.;
+  fV1y = 0.;
+  fV1z = 0.;
+  fV2x = 0.;
+  fV2y = 0.;
+  fV2z = 0.;
+  fDCA = 0.;
+
+  fPx[0] = 0.;
+  fPy[0] = 0.;
+  fPz[0] = 0.;
+  fPx[1] = 0.;
+  fPy[1] = 0.;
+  fPz[1] = 0.;
+
+  fd0[0] = 0.;
+  fd0[1] = 0.;
+
+  fPdg[0] = 0;
+  fPdg[1] = 0;
+  fMum[0] = 0;
+  fMum[1] = 0;
+  fGMum[0] = 0;
+  fGMum[1] = 0;
+
+  fTagPi[0] = 0.;
+  fTagPi[1] = 0.;
+  fTagKa[0] = 0.;
+  fTagKa[1] = 0.;
+  fTagNid[0] = 0.;
+  fTagNid[1] = 0.;
+
+  fWgtJPsi=0;
+
+}
+//----------------------------------------------------------------------------
+AliBtoJPSItoEle::AliBtoJPSItoEle(Int_t ev,Int_t trkNum[2],
+                      Double_t v1[3],Double_t v2[3], 
+                      Double_t dca,
+                      Double_t mom[6],Double_t d0[2]) {
+  // Constructor
+
+  fSignal = kFALSE;
+  fJpsiPrimary = kFALSE;
+
+  fEvent = ev;
+  fTrkNum[0] = trkNum[0];
+  fTrkNum[1] = trkNum[1];
+
+  fV1x = v1[0];
+  fV1y = v1[1];
+  fV1z = v1[2];
+  fV2x = v2[0];
+  fV2y = v2[1];
+  fV2z = v2[2];
+  fDCA = dca;
+
+  fPx[0] = mom[0];
+  fPy[0] = mom[1];
+  fPz[0] = mom[2];
+  fPx[1] = mom[3];
+  fPy[1] = mom[4];
+  fPz[1] = mom[5];
+
+  fd0[0] = d0[0];
+  fd0[1] = d0[1];
+
+  fPdg[0] = 0;
+  fPdg[1] = 0;
+  fMum[0] = 0;
+  fMum[1] = 0;
+  fGMum[0] = 0;
+  fGMum[1] = 0;
+
+  fTagPi[0]  = 0.;
+  fTagPi[1]  = 0.;
+  fTagKa[0]  = 0.;
+  fTagKa[1]  = 0.;
+  fTagNid[0] = 0.;
+  fTagNid[1] = 0.;
+
+  fWgtJPsi=0;
+}
+//----------------------------------------------------------------------------
+AliBtoJPSItoEle::~AliBtoJPSItoEle() {}
+//____________________________________________________________________________
+AliBtoJPSItoEle::AliBtoJPSItoEle( const AliBtoJPSItoEle& btoJpsi):TObject(btoJpsi) {
+  // dummy copy constructor
+}
+//----------------------------------------------------------------------------
+void AliBtoJPSItoEle::ApplyPID(TString pidScheme) {
+  // Applies particle identification
+
+  if(!pidScheme.CompareTo("TRDTPCparam")  && fPdg[0]==0) {
+    printf("AliBtoJPSItoEle::ApplyPID :\n Warning: TRD-TPC parameterized PID can be used only for simulation!\n"); 
+    return;
+  }
+
+  if(!pidScheme.CompareTo("TRDTPCparam")) {
+    // tagging of the positive track
+    if(TMath::Abs(fPdg[0])==11) { // electron
+      fTagEl[0] = 0.81;
+      fTagNid[0] = 1.-fTagEl[0];
+    }
+    else if(TMath::Abs(fPdg[0])==211) { // pion
+      fTagEl[0]   = TRDTPCCombinedPIDParametrization(PChild(0));
+      fTagNid[0] = 1.-fTagEl[0];
+    } 
+    else { // all the others 
+      fTagEl[0]  = 0.;
+      fTagNid[0] = 1.;
+    } 
+    // tagging of the negative track
+    if(TMath::Abs(fPdg[1])==11) { // electron
+      fTagEl[1] = 0.81;
+      fTagNid[1] = 1.-fTagEl[1];
+    }
+    else if(TMath::Abs(fPdg[1])==211) { // pion
+      fTagEl[1]   = TRDTPCCombinedPIDParametrization(PChild(1));
+      fTagNid[1] = 1.-fTagEl[1];
+    }
+    else { // all the others
+      fTagEl[1]  = 0.;
+      fTagNid[1] = 1.;
+    }
+  }
+
+  if(!pidScheme.CompareTo("ESDCombinedPID")) {
+    fTagEl[0]=fPIDrespEl[0];
+    fTagEl[1]=fPIDrespEl[1];
+    fTagNid[0] = 1.-fTagEl[0];
+    fTagNid[1] = 1.-fTagEl[1];
+  }
+  return;
+}
+//----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEle::ChildrenRelAngle() const {
+  // relative angle between K and pi
+
+  TVector3 mom0(fPx[0],fPy[0],fPz[0]);
+  TVector3 mom1(fPx[1],fPy[1],fPz[1]);
+
+  Double_t angle = mom0.Angle(mom1);
+
+  return angle; 
+}
+//----------------------------------------------------------------------------
+void AliBtoJPSItoEle::ComputeWgts() {
+  // calculate the weights for PID
+
+
+  // assignement of the weights from PID
+  fWgtJPsi    = fTagEl[0]*fTagEl[1]; // both assumed to be electrons 
+
+  
+  // if(fWgtJPsi<0.) cerr<<"AliBtoJPSItoEle::ComputeWgts()  Negative weight!!!\n";
+  
+
+  return;
+}
+//----------------------------------------------------------------------------
+void AliBtoJPSItoEle::CorrectWgt4BR(Double_t factor) {
+  // correct weights of background from charm 
+
+  fWgtJPsi    *= factor;
+
+  return;
+}
+//----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEle::CosPointing() const {
+  // cosine of pointing angle in space
+
+  TVector3 mom(Px(),Py(),Pz());
+  TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
+
+  Double_t pta = mom.Angle(flight);
+
+  return TMath::Cos(pta); 
+}
+//----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEle::CosPointingXY() const {
+  // cosine of pointing angle in transverse plane
+
+  TVector3 momXY(Px(),Py(),0.);
+  TVector3 flightXY(fV2x-fV1x,fV2y-fV1y,0.);
+
+  Double_t ptaXY = momXY.Angle(flightXY);
+
+  return TMath::Cos(ptaXY); 
+}
+//----------------------------------------------------------------------------
+void AliBtoJPSItoEle::CosThetaStar(Double_t &ctsJPsi) const {
+  // cosine of decay angle in the J/Psi rest frame (of the negative electron)
+
+  Double_t pStar = TMath::Sqrt(TMath::Power(kMJPsi*kMJPsi-2.*kMe*kMe,2.)-4.*kMe*kMe*kMe*kMe)/(2.*kMJPsi);
+
+  Double_t beta = P()/Energy();
+  Double_t gamma = Energy()/kMJPsi;
+
+  ctsJPsi = (Ql(1)/gamma-beta*TMath::Sqrt(pStar*pStar+kMe*kMe))/pStar;
+  // if(ctsJPsi > 1.)  { cerr<<"AliBtoJPSItoEle::CosThetaStar: > 1 "<<ctsJPsi<<"!\n"; }
+  // if(ctsJPsi < -1.) { cerr<<"AliBtoJPSItoEle::CosThetaStar: < -1 "<<ctsJPsi<<"!\n"; }
+
+  return;
+}
+//----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEle::Eta() const {
+  // pseudorapidity of the J/Psi
+
+  Double_t theta = TMath::Pi()/2.-TMath::ATan2(Pz(),Pt());
+  Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
+  return eta;
+}
+//----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEle::EtaChild(Int_t child) const {
+  // pseudorapidity of the decay tracks
+
+  Double_t theta = TMath::Pi()/2.-TMath::ATan2(fPz[child],PtChild(child));
+  Double_t eta = -TMath::Log(TMath::Tan(theta/2.));
+  return eta;
+}
+//----------------------------------------------------------------------------
+void AliBtoJPSItoEle::GetWgts(Double_t &WgtJPsi) 
+  const {
+  // returns the weights for pid
+
+    WgtJPsi    = fWgtJPsi; 
+
+  return;
+}
+//----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEle::ImpPar() const {
+  // J/Psi impact parameter in the bending plane
+  
+    Double_t k = -(fV2x-fV1x)*Px()-(fV2y-fV1y)*Py();
+    k /= Pt()*Pt();
+    Double_t dx = fV2x-fV1x+k*Px();
+    Double_t dy = fV2y-fV1y+k*Py();
+    Double_t absDD = TMath::Sqrt(dx*dx+dy*dy);
+    TVector3 mom(Px(),Py(),Pz());
+    TVector3 flight(fV2x-fV1x,fV2y-fV1y,fV2z-fV1z);
+    TVector3 cross = mom.Cross(flight);
+    return (cross.Z()>0. ? absDD : -absDD);
+}
+//----------------------------------------------------------------------------
+void AliBtoJPSItoEle::InvMass(Double_t &mJPsi) const {
+  // invariant mass as J/Psi
+
+  Double_t energy[2];
+
+  // J/psi -> e- e+
+  energy[1] = TMath::Sqrt(kMe*kMe+PChild(1)*PChild(1));
+  energy[0] = TMath::Sqrt(kMe*kMe+PChild(0)*PChild(0));
+
+  mJPsi = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-P()*P());
+    
+  return;
+
+}
+//----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEle::Ql(Int_t child) const {
+  // longitudinal momentum of decay tracks w.r.t. to J/Psi momentum
+
+  Double_t qL;
+  TVector3 mom(fPx[child],fPy[child],fPz[child]);
+  TVector3 momJPsi(Px(),Py(),Pz());
+
+  qL = mom.Dot(momJPsi)/momJPsi.Mag();
+
+  return qL ;
+}
+//----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEle::Qt() const {
+  // transverse momentum of decay tracks w.r.t. to JPsi momentum  
+
+  TVector3 mom0(fPx[0],fPy[0],fPz[0]);
+  TVector3 momJPsi(Px(),Py(),Pz());
+
+  return mom0.Perp(momJPsi);
+}
+//----------------------------------------------------------------------------
+Bool_t AliBtoJPSItoEle::Select(const Double_t* cuts,Int_t& okB) 
+  const {
+//
+// This function compares the B candidates with a set of cuts:
+//
+// cuts[0] = inv. mass half width [GeV]   
+// cuts[1] = dca [micron]
+// cuts[2] = cosThetaStar 
+// cuts[3] = pTP [GeV/c]
+// cuts[4] = pTN [GeV/c]
+// cuts[5] = d0P [micron]   upper limit!
+// cuts[6] = d0N [micron]  upper limit!
+// cuts[7] = d0d0 [micron^2]
+// cuts[8] = cosThetaPoint
+//
+// If the candidate doesn't pass the cuts it sets the weight to 0
+// and return kFALSE
+//
+  Double_t mJPsi,ctsJPsi;
+  okB=1; 
+
+  if(PtChild(1) < cuts[3] || PtChild(0) < cuts[4]) okB = 0;
+  if(!okB) return kFALSE;
+
+  if(TMath::Abs(Getd0Child(1)) > cuts[5] || 
+     TMath::Abs(Getd0Child(0)) > cuts[6]) okB = 0;
+  if(!okB) return kFALSE;
+
+  if(GetDCA() > cuts[1]) { okB = 0; return kFALSE; }
+
+  InvMass(mJPsi);
+  if(TMath::Abs(mJPsi-kMJPsi)    > cuts[0]) okB = 0;
+  if(!okB) return kFALSE;
+
+  CosThetaStar(ctsJPsi);
+  if(TMath::Abs(ctsJPsi)    > cuts[2]) okB = 0;
+  if(!okB) return kFALSE;
+
+  if(ProdImpParams() > cuts[7]) { okB = 0; return kFALSE; }
+
+  if(CosPointing()   < cuts[8]) { okB = 0; return kFALSE; }
+
+  return kTRUE;
+}
+//-----------------------------------------------------------------------------
+void AliBtoJPSItoEle::SetPIDresponse(Double_t resp0[5],Double_t resp1[5]) {
+  // Set combined PID detector response probabilities
+
+  fPIDrespEl[0] = resp0[0];
+  fPIDrespEl[1] = resp1[0];
+  fPIDrespMu[0] = resp0[1];
+  fPIDrespMu[1] = resp1[1];
+  fPIDrespPi[0] = resp0[2];
+  fPIDrespPi[1] = resp1[2];
+  fPIDrespKa[0] = resp0[3];
+  fPIDrespKa[1] = resp1[3];
+  fPIDrespPr[0] = resp0[4];
+  fPIDrespPr[1] = resp1[4];
+
+  return;
+} 
+//-----------------------------------------------------------------------------
+void AliBtoJPSItoEle::GetPIDresponse(Double_t resp0[5],Double_t resp1[5]) const {
+  // Get combined PID detector response probabilities
+                                                                                                                          
+  resp0[0] = fPIDrespEl[0];
+  resp1[0] = fPIDrespEl[1]; 
+  resp0[1] = fPIDrespMu[0];
+  resp1[1] = fPIDrespMu[1];
+  resp0[2] = fPIDrespPi[0];
+  resp1[2] = fPIDrespPi[1];
+  resp0[3] = fPIDrespKa[0];
+  resp1[3] = fPIDrespKa[1];
+  resp0[4] = fPIDrespPr[0];
+  resp1[4] = fPIDrespPr[1];
+                                                                                                                          
+  return;
+}
+//----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEle::TRDTPCCombinedPIDParametrization(Double_t p) const {
+
+
+  // a first raw parametrization of the probability to misidentify a charged pion as electron as a 
+  // function of the momentum, as given by the combined TPC and TRD response. 
+  //  PID cuts are set such that the probability for correct electron id is 90% in each of the two 
+  //    detectors
+  
+// first estimate based on parameterization in the B-> single electron analysis
+  Double_t value=0;
+  Double_t p1 =11.;
+  Double_t p2=0.00007;
+  Double_t p3=0.007;
+  value=p2+p3*(1.-exp(-TMath::Power(p/p1,4.)));
+
+// Better estimation based on TRD test beam (as presented by Andrea at Munster)
+  value/=0.01; // here remove from TPC+TRD the TRD contribution estimated to be 0.01
+  if (p<10.) value*=(1.32-0.18*p+0.076*p*p-0.0037*p*p*p)/100.;
+  if (p>10.) value*=(0.48+0.287*p)/100.;
+
+  return value;
+}
+//----------------------------------------------------------------------------
+
+
+
+
+
diff --git a/PWG3/AliBtoJPSItoEle.h b/PWG3/AliBtoJPSItoEle.h
new file mode 100644 (file)
index 0000000..9cc483b
--- /dev/null
@@ -0,0 +1,151 @@
+#ifndef AliBtoJPSItoEle_H
+#define AliBtoJPSItoEle_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                          Class AliBtoJPSItoEle
+//                 Reconstructed B -> J\PSI+X --> e+ e- *X  class
+//      
+// Note: the two decay tracks are labelled: 0 (positive electron)
+//                                          1 (negative electron)
+//
+//         Origin: G.E. Bruno    giuseppe.bruno@ba.infn.it                  
+//          based on Class for charm golden channel (D0->Kpi) (thanks to Andrea Dainese!)
+//-------------------------------------------------------------------------
+
+#include <TObject.h>
+#include <TMath.h>
+
+//----------------------------------------------------------------------------
+//     Some constants (masses)
+//
+// particle masses
+const Double_t kMJPsi = 3.096916;     // J/Psi mass
+const Double_t kMe    = 0.000510998902;  // elec  mass
+
+
+
+
+//-----------------------------------------------------------------------------
+class AliBtoJPSItoEle : public TObject {
+ public:
+  //
+  AliBtoJPSItoEle();
+  AliBtoJPSItoEle(Int_t ev,Int_t trkNum[2],  
+            Double_t v1[3],Double_t v2[3],Double_t dca,
+            Double_t mom[6],Double_t d0[2]);
+  virtual ~AliBtoJPSItoEle();
+  AliBtoJPSItoEle(const AliBtoJPSItoEle& btoJpsi);
+
+  Double_t Alpha() const { return (Ql(0)-Ql(1))/(Ql(0)+Ql(1)); }
+  void     ApplyPID(TString pidScheme="TRDTPCparam");
+  Double_t ChildrenRelAngle() const; 
+  void     ComputeWgts();
+  void     CorrectWgt4BR(Double_t factor);
+  Double_t CosPointing() const;
+  Double_t CosPointingXY() const;
+  void     CosThetaStar(Double_t &ctsJPSI) const;
+  Double_t Ct() const {return Length()*kMJPsi/P();}
+  Double_t Energy() const { return TMath::Sqrt(P()*P()+kMJPsi*kMJPsi); }
+  Double_t Eta() const;
+  Double_t EtaChild(Int_t child) const;
+  Int_t    EventNo() const {return TMath::Abs(fEvent);}
+  Double_t GetDCA() const { return 10000.*fDCA; }
+  Int_t    GetTrkNum(Int_t child) const { return fTrkNum[child]; }
+  Double_t Getd0Child(Int_t child) const { return fd0[child]; }
+  Int_t    GetPdgChild(Int_t child) const { return fPdg[child]; }
+  Int_t    GetPdgMum(Int_t child) const {return fMum[child]; }
+  Int_t    GetPdgGMum(Int_t child) const {return fGMum[child]; }
+  void     GetPIDresponse(Double_t resp0[5],Double_t resp1[5]) const; 
+  void     GetWgts(Double_t &WgtJPsi) const;
+  void     GetPrimaryVtx(Double_t vtx[3]) const 
+    { vtx[0]=fV1x; vtx[1]=fV1y; vtx[2]=fV1z; }
+  void     GetSecondaryVtx(Double_t vtx[3]) const 
+    { vtx[0]=fV2x; vtx[1]=fV2y; vtx[2]=fV2z; }
+
+  Double_t ImpPar() const;
+  void     InvMass(Double_t &mJPsi) const;
+  Bool_t   IsSignal() const { if(fSignal) return kTRUE; return kFALSE; } 
+  Bool_t   IsJpsiPrimary() const { if(fJpsiPrimary) return kTRUE; return kFALSE; } 
+  Double_t Length() const
+    { return TMath::Sqrt((fV1x-fV2x)*(fV1x-fV2x)
+                        +(fV1y-fV2y)*(fV1y-fV2y)+(fV1z-fV2z)*(fV1z-fV2z)); }
+  Double_t P()  const { return TMath::Sqrt(Pt()*Pt()+Pz()*Pz()); } 
+  Double_t PChild(Int_t child)   const { return TMath::Sqrt(fPx[child]*fPx[child]+fPy[child]*fPy[child]+fPz[child]*fPz[child]); } 
+  Double_t ProdImpParams() const { return fd0[0]*fd0[1]; } 
+  Double_t Pt() const { return TMath::Sqrt(Px()*Px()+Py()*Py()); }
+  Double_t PtChild(Int_t child)  const { return TMath::Sqrt(fPx[child]*fPx[child]+fPy[child]*fPy[child]); } 
+  Double_t Px() const { return (fPx[0]+fPx[1]); }
+  Double_t Py() const { return (fPy[0]+fPy[1]); }
+  Double_t Pz() const { return (fPz[0]+fPz[1]); }
+  Double_t Ql(Int_t child) const;
+  Double_t Qt() const;
+  Double_t Rapidity() const { return 0.5*TMath::Log((Energy()+Pz())/(Energy()-Pz()+1.e-13)); }
+  Bool_t   Select(const Double_t* cuts,Int_t& okB) const;
+  void     SetPrimaryVtx(Double_t vtx[3]) 
+    { fV1x=vtx[0]; fV1y=vtx[1]; fV1z=vtx[2]; }
+  void     SetSignal() { fSignal =  kTRUE; }
+  void     SetJpsiPrimary() { fJpsiPrimary =  kTRUE; }
+  void     SetTOFmasses(Double_t mass[2]) 
+    { fTOFmass[0]=mass[0]; fTOFmass[1]=mass[1]; }
+  void     SetPIDresponse(Double_t resp0[5],Double_t resp1[5]); 
+  void     SetPdgCodes(Int_t pdg[2]) {fPdg[0]=pdg[0];fPdg[1]=pdg[1]; }
+  void     SetMumPdgCodes(Int_t mum[2]) {fMum[0]=mum[0];fMum[1]=mum[1]; }
+  void     SetGMumPdgCodes(Int_t gmum[2]) {fGMum[0]=gmum[0];fGMum[1]=gmum[1]; }
+  Double_t TRDTPCCombinedPIDParametrization(Double_t p) const;
+  //
+ private:
+  //
+  Bool_t   fSignal; // TRUE if signal, FALSE if background (for simulation) 
+                    // (background are both combinatorial and primary J/psi)
+  Bool_t   fJpsiPrimary; // TRUE if the current candidate is a primary J/psi, FALSE otherway (for simulation) 
+  Int_t    fEvent;  // number of the event this B comes from
+                    // -1 if the B comes from ev. mixing
+
+  Int_t fTrkNum[2]; // numbers of the two decay tracks
+
+  Double_t fV1x; // X-position of the primary vertex of the event
+  Double_t fV1y; // Y-position of the primary vertex of the event
+  Double_t fV1z; // Z-position of the primary vertex of the event
+  Double_t fV2x; // X-position of the reconstructed secondary vertex
+  Double_t fV2y; // Y-position of the reconstructed secondary vertex
+  Double_t fV2z; // Z-position of the reconstructed secondary vertex
+  Double_t fDCA; // DCA of the two tracks
+
+  Double_t fPx[2];  // X,Y,Z
+  Double_t fPy[2];  // momenta of the two tracks
+  Double_t fPz[2];  // at the reconstructed vertex  
+
+  Double_t fd0[2];  //  impact parameters in the bending plane
+
+  Int_t fPdg[2];  // PDG codes of the two tracks       (for sim.)
+  Int_t fMum[2];  // PDG codes of the mothers          (for sim.)
+  Int_t fGMum[2]; // PDG codes of the grand-mothers    (for sim.)
+
+  Double_t fTagEl[2];  // probability to be tagged as electron 
+  Double_t fTagPi[2];  // probability to be tagged as pion 
+  Double_t fTagKa[2];  // probability to be tagged as kaon 
+  Double_t fTagPr[2];  // probability to be tagged as proton 
+  Double_t fTagNid[2]; // probability to be tagged as "non-identified" 
+
+  Double_t fPIDrespEl[2]; // det. response to be electron
+  Double_t fPIDrespMu[2]; // det. response to be muon
+  Double_t fPIDrespPi[2]; // det. response to be pion
+  Double_t fPIDrespKa[2]; // det. response to be kaon
+  Double_t fPIDrespPr[2]; // det. response to be proton
+  Double_t fTOFmass[2]; // mass estimated by the TOF (-1000. if track not reached TOF)
+  Double_t fWgtJPsi; // weights for the pair
+
+  ClassDef(AliBtoJPSItoEle,1)  // Reconstructed B candidate class
+};
+
+#endif
+
+
+
+
+
+
+
+
diff --git a/PWG3/AliBtoJPSItoEleAnalysis.cxx b/PWG3/AliBtoJPSItoEleAnalysis.cxx
new file mode 100644 (file)
index 0000000..c9c2a54
--- /dev/null
@@ -0,0 +1,721 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//----------------------------------------------------------------------------
+//    Implementation of the BtoJPSItoEle reconstruction and analysis class
+// Note: the two decay tracks are labelled: 0 (positive track)
+//                                          1 (negative track)
+// An example of usage can be found in the macro AliBtoJPSItoEleTest.C
+//            Origin: G.E. Bruno    giuseppe.bruno@ba.infn.it            
+//  based on Class for charm golden channel (D0->Kpi)
+//----------------------------------------------------------------------------
+#include <TKey.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TString.h>
+#include <TSystem.h>
+#include <TParticle.h>
+#include "AliESD.h"
+#include "AliMC.h"
+#include "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliVertexerTracks.h"
+#include "AliESDVertex.h"
+#include "AliESDv0.h"
+#include "AliBtoJPSItoEle.h"
+#include "AliBtoJPSItoEleAnalysis.h"
+#include "AliLog.h"
+
+typedef struct {
+  Int_t lab;
+  Int_t pdg;
+  Int_t mumlab;
+  Int_t mumpdg;
+  Int_t gmumlab;
+  Int_t gmumpdg;
+  Int_t mumprongs;
+  Float_t Vx,Vy,Vz;
+  Float_t Px,Py,Pz;
+} REFTRACK;
+
+ClassImp(AliBtoJPSItoEleAnalysis)
+
+//----------------------------------------------------------------------------
+AliBtoJPSItoEleAnalysis::AliBtoJPSItoEleAnalysis() {
+  // Default constructor
+
+  SetPtCut();
+  Setd0Cut();
+  SetPidCut();
+  SetMassCut();
+  SetBCuts();
+  SetVertex1();
+  SetPID();
+  fVertexOnTheFly = kFALSE;
+  fSim = kFALSE;
+  fOnlySignal = kFALSE;
+  fOnlyPrimaryJpsi = kFALSE;
+}
+//----------------------------------------------------------------------------
+AliBtoJPSItoEleAnalysis::~AliBtoJPSItoEleAnalysis() {}
+//----------------------------------------------------------------------------
+void AliBtoJPSItoEleAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
+  // select candidates that pass fBCuts and write them to a new file
+
+  TFile *inFile = TFile::Open(inName);
+
+  TTree *treeBin=(TTree*)inFile->Get("TreeB");
+  AliBtoJPSItoEleAnalysis *inAnalysis = (AliBtoJPSItoEleAnalysis*)inFile->Get("BtoJPSItoEleAnalysis");
+  printf("+++\n+++  I N P U T   S T A T U S:\n+++\n");
+  inAnalysis->PrintStatus();
+
+
+  AliBtoJPSItoEle *d = 0; 
+  treeBin->SetBranchAddress("BtoJPSItoEle",&d);
+  Int_t entries = (Int_t)treeBin->GetEntries();
+
+  printf("+++\n+++ Number of B candidates in input tree:  %d\n+++\n",entries);
+
+  TTree *treeBout = new TTree("TreeB","Tree with selected B candidates");
+  treeBout->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&d,200000,0);
+
+
+  Int_t okB=0;
+  Int_t nSel = 0;
+
+  for(Int_t i=0; i<entries; i++) {
+    // get event from tree
+    treeBin->GetEvent(i);
+
+    if(fSim && fOnlySignal && !d->IsSignal()) continue; 
+
+    // check if candidate passes selection (as B or Bbar)
+    if(d->Select(fBCuts,okB)) {
+      nSel++;
+      treeBout->Fill();
+    }
+
+  }
+
+  AliBtoJPSItoEleAnalysis *outAnalysis = (AliBtoJPSItoEleAnalysis*)inAnalysis->Clone("BtoJPSItoEleAnalysis");
+  outAnalysis->SetBCuts(fBCuts);
+  printf("------------------------------------------\n");
+  printf("+++\n+++  O U T P U T   S T A T U S:\n+++\n");
+  outAnalysis->PrintStatus();
+
+  printf("+++\n+++ Number of B mesons in output tree:  %d\n+++\n",nSel);
+
+  TFile* outFile = new TFile(outName,"recreate");
+  treeBout->Write();
+  outAnalysis->Write();
+  outFile->Close();
+
+  return;
+}
+//----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEleAnalysis::CalculateTOFmass(Double_t mom,Double_t length,
+                                             Double_t time) const {
+  // calculated the mass from momentum, track length from vertex to TOF
+  // and time measured by the TOF
+  if(length==0.) return -1000.;
+  Double_t a = time*time/length/length;
+  if(a > 1.) {
+    a = TMath::Sqrt(a-1.);
+  } else {
+    a = -TMath::Sqrt(1.-a);
+  }
+
+  return mom*a;
+}
+//----------------------------------------------------------------------------
+void AliBtoJPSItoEleAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
+                                       const Char_t *outName) {
+  // Find candidates and calculate parameters
+
+
+  TString esdName="AliESDs.root";
+  if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
+    printf("AliBtoJPSItoEleAnalysis::FindCandidatesESD(): No ESDs file found!\n"); 
+    return;
+  }
+
+  TString outName1=outName;
+  TString outName2="nTotEvents.dat";
+
+  Int_t    nTotEv=0,nBrec=0,nBrec1ev=0;
+  Double_t dca;
+  Double_t v2[3],mom[6],d0[2];
+  Int_t    iTrkP,iTrkN,trkEntries;
+  Int_t    nTrksP=0,nTrksN=0;
+  Int_t    trkNum[2];
+  Double_t tofmass[2];
+  Int_t    okB=0;
+  AliESDtrack *postrack = 0;
+  AliESDtrack *negtrack = 0;
+
+  // create the AliVertexerTracks object
+  // (it will be used only if fVertexOnTheFly=kTrue)
+  AliVertexerTracks *vertexer1 = new AliVertexerTracks;
+  if(fVertexOnTheFly) {
+    // open the mean vertex
+    TFile *invtx = new TFile("AliESDVertexMean.root");
+    AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
+    invtx->Close();
+    vertexer1->SetVtxStart(initVertex);
+    delete invtx;
+  }
+  Int_t  skipped[2];
+  Bool_t goodVtx1;
+  
+  // create tree for reconstructed decayes
+  AliBtoJPSItoEle *ioBtoJPSItoEle=0;
+  TTree *treeB = new TTree("TreeB","Tree with candidates");
+  treeB->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&ioBtoJPSItoEle,200000,0);
+
+  // open file with tracks
+  TFile *esdFile = TFile::Open(esdName.Data());
+  AliESD* event = new AliESD;
+  TTree* tree = (TTree*) esdFile->Get("esdTree");
+  if(!tree) {
+    Error("FindCandidatesESD", "no ESD tree found");
+    return;
+  }
+  tree->SetBranchAddress("ESD",&event);
+
+  // loop on events in file
+  for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
+    if(iEvent > evLast) break;
+    tree->GetEvent(iEvent);
+    Int_t ev = (Int_t)event->GetEventNumberInFile();
+    printf("--- Finding B -> JPSI -> e+ e- in event %d\n",ev);
+    // count the total number of events
+    nTotEv++;
+
+    trkEntries = (Int_t)event->GetNumberOfTracks();
+    printf(" Number of tracks: %d\n",trkEntries);
+    if(trkEntries<2) continue;
+
+    // retrieve primary vertex from file
+    if(!event->GetPrimaryVertex()) { 
+      printf(" No vertex\n");
+      continue;
+    }
+    event->GetPrimaryVertex()->GetXYZ(fV1);
+
+    // call function which applies sigle-track selection and
+    // separetes positives and negatives
+    TObjArray trksP(trkEntries/2);
+    Int_t    *trkEntryP   = new Int_t[trkEntries];
+    TObjArray trksN(trkEntries/2);
+    Int_t    *trkEntryN   = new Int_t[trkEntries];
+    TTree *trkTree = new TTree();
+    SelectTracks(event,trksP,trkEntryP,nTrksP,
+                      trksN,trkEntryN,nTrksN);
+
+    printf(" pos. tracks: %d    neg .tracks: %d\n",nTrksP,nTrksN);
+
+
+    nBrec1ev = 0;
+
+    // LOOP ON  POSITIVE  TRACKS
+    for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
+      if(iTrkP%1==0) printf("  Processing positive track number %d of %d\n",iTrkP,nTrksP);
+         
+      // get track from track array
+      postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
+      trkNum[0] = trkEntryP[iTrkP];      
+
+      // LOOP ON  NEGATIVE  TRACKS
+      for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
+
+       // get track from tracks array
+       negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
+       trkNum[1] = trkEntryN[iTrkN];      
+
+       //
+       // ----------- DCA MINIMIZATION ------------------
+       //
+       // find the DCA and propagate the tracks to the DCA
+       Double_t b=event->GetMagneticField(); 
+       AliESDtrack nt(*negtrack), pt(*postrack);
+       dca = nt.PropagateToDCA(&pt,b);
+
+       // define the AliESDv0 object
+       AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
+         
+       // get position of the secondary vertex
+       vertex2.GetXYZ(v2[0],v2[1],v2[2]);
+        vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]);
+        vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]);
+       // impact parameters of the tracks w.r.t. the primary vertex
+       // d0[0] =  10000.*pt.GetD(fV1[0],fV1[1],b);
+       // d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
+       goodVtx1 = kTRUE;
+       
+       // no vertexing if DeltaMass > fMassCut 
+       if(fVertexOnTheFly) {
+         goodVtx1 = kFALSE;
+         if(SelectInvMass(mom)) {
+           // primary vertex from *other* tracks in the event
+           skipped[0] = trkEntryP[iTrkP];
+           skipped[1] = trkEntryN[iTrkN];
+           vertexer1->SetSkipTracks(2,skipped);
+           AliESDVertex *vertex1onfly = 
+             (AliESDVertex*)vertexer1->FindPrimaryVertex(event); 
+           if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
+           vertex1onfly->GetXYZ(fV1);
+           //vertex1onfly->PrintStatus();
+           delete vertex1onfly;        
+         }
+       }
+       
+       // impact parameters of the tracks w.r.t. the primary vertex
+       d0[0] =  10000.*pt.GetD(fV1[0],fV1[1],b);
+       d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
+
+       // create the object AliBtoJPSItoEle
+       AliBtoJPSItoEle theB(ev,trkNum,fV1,v2,dca,mom,d0);
+       // select B's
+       if(goodVtx1 && theB.Select(fBCuts,okB)) {
+         // get PID info from ESD
+         AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
+         Double_t esdpid0[5];
+         t0->GetESDpid(esdpid0);
+         if(t0->GetStatus()&AliESDtrack::kTOFpid) {
+           tofmass[0] = CalculateTOFmass(t0->GetP(),
+                                         t0->GetIntegratedLength(),
+                                         t0->GetTOFsignal());
+         } else {
+           tofmass[0] = -1000.;
+         }
+         AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
+         Double_t esdpid1[5];
+         t1->GetESDpid(esdpid1);
+         if(t1->GetStatus()&AliESDtrack::kTOFpid) {
+           tofmass[1] = CalculateTOFmass(t1->GetP(),
+                                         t1->GetIntegratedLength(),
+                                         t1->GetTOFsignal());
+         } else {
+           tofmass[1] = -1000.;
+         }
+
+         theB.SetPIDresponse(esdpid0,esdpid1);
+         theB.SetTOFmasses(tofmass);
+
+         // fill the tree
+         ioBtoJPSItoEle=&theB;
+         treeB->Fill();
+
+         nBrec++; nBrec1ev++;
+         ioBtoJPSItoEle=0; 
+       }
+       
+       negtrack = 0;
+      } // loop on negative tracks
+      postrack = 0;
+    }   // loop on positive tracks
+    
+    delete [] trkEntryP;
+    delete [] trkEntryN;
+    delete trkTree;
+
+    printf(" Number of B candidates: %d\n",nBrec1ev);
+  }    // loop on events in file
+
+
+  printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
+  printf("\n+++\n+++ Total number of B candidates: %d\n+++\n",nBrec);
+
+  delete vertexer1;
+
+  esdFile->Close();
+
+  // create a copy of this class to be written to output file
+  AliBtoJPSItoEleAnalysis *copy = (AliBtoJPSItoEleAnalysis*)this->Clone("BtoJPSItoEleAnalysis");
+
+  // add PDG codes to decay tracks in found candidates (in simulation mode)
+  // and store tree in the output file
+  if(!fSim) {
+    TFile *outroot = new TFile(outName1.Data(),"recreate");
+    treeB->Write();
+    copy->Write();
+    outroot->Close();
+    delete outroot;
+  } else {
+    printf(" Now adding information from simulation (PDG codes) ...\n");
+    TTree *treeBsim = new TTree("TreeB","Tree with B candidates");
+    SimulationInfo(treeB,treeBsim);
+    delete treeB;
+    TFile *outroot = new TFile(outName1.Data(),"recreate");
+    treeBsim->Write();
+    copy->Write();
+    outroot->Close();
+    delete outroot;
+  }
+
+  // write to a file the total number of events
+  FILE *outdat = fopen(outName2.Data(),"w");
+  fprintf(outdat,"%d\n",nTotEv);
+  fclose(outdat);
+
+  return;
+}
+//-----------------------------------------------------------------------------
+void AliBtoJPSItoEleAnalysis::PrintStatus() const {
+  // Print parameters being used
+
+  printf("Preselections:\n");
+  printf("    fPtCut   = %f GeV\n",fPtCut);
+  printf("    fd0Cut   = %f micron\n",fd0Cut);
+  printf("    fMassCut = %f GeV\n",fMassCut);
+  printf("    fPidCut  > %f \n",fPidCut);
+  if(fVertexOnTheFly) printf("Primary vertex on the fly\n");
+  if(fSim) { 
+    printf("Simulation mode\n");
+    if(fOnlySignal && !(fOnlyPrimaryJpsi)) printf("  Only signal goes to file\n");
+    if(fOnlyPrimaryJpsi && !(fOnlySignal)) printf("  Only primary Jpsi go to file\n");
+    if(fOnlyPrimaryJpsi && fOnlySignal) printf("  Both signal and primary Jpsi go to file\n");
+  }
+  printf("Cuts on candidates:\n");
+  printf("    |M-MJPsi| [GeV]  < %f\n",fBCuts[0]);
+  printf("    dca    [micron]  < %f\n",fBCuts[1]);
+  printf("    cosThetaStar     < %f\n",fBCuts[2]);
+  printf("    pTP     [GeV]    > %f\n",fBCuts[3]);
+  printf("    pTN     [GeV]    > %f\n",fBCuts[4]);
+  printf("    |d0P|  [micron]  < %f\n",fBCuts[5]);
+  printf("    |d0N|  [micron]  < %f\n",fBCuts[6]);
+  printf("    d0d0  [micron^2] < %f\n",fBCuts[7]);
+  printf("    cosThetaPoint    > %f\n",fBCuts[8]);
+
+  return;
+}
+//-----------------------------------------------------------------------------
+Bool_t AliBtoJPSItoEleAnalysis::SelectInvMass(const Double_t p[6]) const {
+  // Apply preselection in the invariant mass of the pair
+
+  Double_t mJPsi = 3.096916;
+  Double_t mel   = 0.00510998902;
+
+  Double_t energy[2];
+  Double_t mom2[2],momTot2;
+
+  mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
+  mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5];
+
+  momTot2 = (p[0]+p[3])*(p[0]+p[3])+
+            (p[1]+p[4])*(p[1]+p[4])+
+            (p[2]+p[5])*(p[2]+p[5]);
+
+  // J/Psi -> e+ e-
+  energy[1] = TMath::Sqrt(mel*mel+mom2[1]);
+  energy[0] = TMath::Sqrt(mel*mel+mom2[0]);
+
+  Double_t minvJPsi = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
+
+  if(TMath::Abs(minvJPsi-mJPsi)  < fMassCut) return kTRUE;
+  return kFALSE;
+}
+//-----------------------------------------------------------------------------
+void AliBtoJPSItoEleAnalysis::SelectTracks(AliESD *event,
+        TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
+        TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
+  // Create two TObjArrays with positive and negative tracks and 
+  // apply single-track preselection
+
+  nTrksP=0,nTrksN=0;
+
+  Int_t entr = event->GetNumberOfTracks();
+  // transfer ITS tracks from ESD to arrays and to a tree
+  for(Int_t i=0; i<entr; i++) {
+
+    AliESDtrack *esdtrack = event->GetTrack(i);
+    UInt_t status = esdtrack->GetStatus();
+
+    if(!(status&AliESDtrack::kITSin)) continue;
+
+    // single track selection
+    if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
+
+    if(esdtrack->GetSign()<0) { // negative track
+      trksN.AddLast(esdtrack);
+      trkEntryN[nTrksN] = i;
+      nTrksN++;
+    } else {                 // positive track
+      trksP.AddLast(esdtrack);
+      trkEntryP[nTrksP] = i;
+      nTrksP++;
+    }
+
+  } // loop on ESD tracks
+
+  return;
+}
+//-----------------------------------------------------------------------------
+void AliBtoJPSItoEleAnalysis::SetBCuts(Double_t cut0,Double_t cut1,
+                                  Double_t cut2,Double_t cut3,Double_t cut4,
+                                  Double_t cut5,Double_t cut6,
+                                  Double_t cut7,Double_t cut8) {
+  // Set the cuts for B selection
+  fBCuts[0] = cut0;
+  fBCuts[1] = cut1;
+  fBCuts[2] = cut2;
+  fBCuts[3] = cut3;
+  fBCuts[4] = cut4;
+  fBCuts[5] = cut5;
+  fBCuts[6] = cut6;
+  fBCuts[7] = cut7;
+  fBCuts[8] = cut8;
+
+  return;
+}
+//-----------------------------------------------------------------------------
+void AliBtoJPSItoEleAnalysis::SetBCuts(const Double_t cuts[9]) {
+  // Set the cuts for B selection
+
+  for(Int_t i=0; i<9; i++) fBCuts[i] = cuts[i];
+
+  return;
+}
+//-----------------------------------------------------------------------------
+Bool_t 
+AliBtoJPSItoEleAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
+  // Check if track passes some kinematical cuts  
+  // Magnetic field "b" (kG)
+
+  if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut) 
+    return kFALSE;
+  if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut) 
+    return kFALSE;
+  //select only tracks with the "combined PID"
+  UInt_t status = trk.GetStatus();
+  if ((status&AliESDtrack::kESDpid)==0) return kTRUE;
+  Double_t r[5];
+  trk.GetESDpid(r);
+  if(r[0] < fPidCut) return kFALSE;
+
+  return kTRUE;
+}
+//----------------------------------------------------------------------------
+void AliBtoJPSItoEleAnalysis::MakeTracksRefFile(AliRun *gAlice,
+                                          Int_t evFirst,Int_t evLast) const {
+  // Create a file with simulation info for the reconstructed tracks
+  
+  TFile *outFile = TFile::Open("BTracksRefFile.root","recreate");
+  TFile *esdFile = TFile::Open("AliESDs.root");
+
+  AliMC *mc = gAlice->GetMCApp();
+  
+  Int_t      label;
+  TParticle *part;  
+  TParticle *mumpart;
+  TParticle *gmumpart;
+  REFTRACK   reftrk;
+  
+  AliESD* event = new AliESD;
+  TTree* tree = (TTree*) esdFile->Get("esdTree");
+  tree->SetBranchAddress("ESD",&event);
+  // loop on events in file
+  for(Int_t iEvent=evFirst; iEvent<tree->GetEntries(); iEvent++) {
+    if(iEvent>evLast) break;
+    tree->GetEvent(iEvent);
+    Int_t ev = (Int_t)event->GetEventNumberInFile();
+
+    gAlice->GetEvent(ev);
+
+    Int_t nentr=(Int_t)event->GetNumberOfTracks();
+
+    // Tree for true track parameters
+    char ttname[100];
+    sprintf(ttname,"Tree_Ref_%d",ev);
+    TTree *reftree = new TTree(ttname,"Tree with true track params");
+    reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:gmumlab:gmumpdg:mumprongs:Vx/F:Vy:Vz:Px:Py:Pz");
+//    reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
+
+    for(Int_t i=0; i<nentr; i++) {
+      AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
+      label = TMath::Abs(esdtrack->GetLabel());
+
+      part = (TParticle*)mc->Particle(label); 
+      reftrk.lab = label;
+      reftrk.pdg = part->GetPdgCode();
+      reftrk.mumlab = part->GetFirstMother();
+      if(part->GetFirstMother()>=0) {
+       mumpart = (TParticle*)gAlice->GetMCApp()->Particle(part->GetFirstMother());
+       reftrk.mumpdg = mumpart->GetPdgCode();
+       reftrk.mumprongs = mumpart->GetNDaughters();
+        reftrk.gmumlab = mumpart->GetFirstMother();
+        if(mumpart->GetFirstMother()>=0) {
+          gmumpart = (TParticle*)gAlice->GetMCApp()->Particle(mumpart->GetFirstMother());
+          reftrk.gmumpdg = gmumpart->GetPdgCode();
+        }
+      } else {
+       reftrk.mumpdg=-1;
+       reftrk.mumprongs=-1;
+       reftrk.gmumpdg=-1;
+        reftrk.gmumlab=part->GetFirstMother(); // If it hasn't any mother, then it has neither Gmother!
+        // reftrk.gmumlab=-1; // If it hasn't any mother, then it has neither Gmother!
+      }
+      reftrk.Vx = part->Vx();
+      reftrk.Vy = part->Vy();
+      reftrk.Vz = part->Vz();
+      reftrk.Px = part->Px();
+      reftrk.Py = part->Py();
+      reftrk.Pz = part->Pz();
+      
+      reftree->Fill();
+
+    } // loop on tracks   
+
+    outFile->cd();
+    reftree->Write();
+
+    delete reftree;
+  } // loop on events
+
+  esdFile->Close();
+  outFile->Close();
+
+  return;
+}
+//-----------------------------------------------------------------------------
+void AliBtoJPSItoEleAnalysis::SimulationInfo(TTree *treeBin,TTree *treeBout) const {
+  // add pdg codes to candidate decay tracks (for sim)
+
+  TString refFileName("BTracksRefFile.root");
+  if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) { 
+    printf("AliBtoJPSItoEleAnalysis::SimulationInfo: no reference file found!\n"); 
+    return;
+  }
+  TFile *refFile = TFile::Open(refFileName.Data());
+
+  Char_t refTreeName[100];
+  Int_t  event;
+  Int_t  pdg[2],mumpdg[2],mumlab[2],gmumpdg[2],gmumlab[2];
+  REFTRACK reftrk;
+
+  // read-in reference tree for event 0 (the only event for Pb-Pb)
+  sprintf(refTreeName,"Tree_Ref_%d",0);
+  TTree *refTree0 = (TTree*)refFile->Get(refTreeName);
+  refTree0->SetBranchAddress("rectracks",&reftrk);
+
+  AliBtoJPSItoEle *theB = 0; 
+  treeBin->SetBranchAddress("BtoJPSItoEle",&theB);
+  treeBout->Branch("BtoJPSItoEle","AliBtoJPSItoEle",&theB,200000,0);
+
+  Int_t entries = (Int_t)treeBin->GetEntries();
+
+  for(Int_t i=0; i<entries; i++) {
+    if(i%100==0) printf("  done %d candidates of %d\n",i,entries);    
+
+    treeBin->GetEvent(i);
+    event = theB->EventNo();
+
+    if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
+      refTree0->GetEvent(theB->GetTrkNum(0));
+      pdg[0]    = reftrk.pdg;
+      mumpdg[0] = reftrk.mumpdg;
+      mumlab[0] = reftrk.mumlab;
+      gmumpdg[0] = reftrk.gmumpdg;
+      gmumlab[0] = reftrk.gmumlab;
+      refTree0->GetEvent(theB->GetTrkNum(1));
+      pdg[1]    = reftrk.pdg;
+      mumpdg[1] = reftrk.mumpdg;
+      mumlab[1] = reftrk.mumlab;
+      gmumpdg[1] = reftrk.gmumpdg;
+      gmumlab[1] = reftrk.gmumlab;
+    } else {
+      sprintf(refTreeName,"Tree_Ref_%d",event);
+      TTree *refTree = (TTree*)refFile->Get(refTreeName);
+      refTree->SetBranchAddress("rectracks",&reftrk);
+      refTree->GetEvent(theB->GetTrkNum(0));
+      pdg[0]    = reftrk.pdg;
+      mumpdg[0] = reftrk.mumpdg;
+      mumlab[0] = reftrk.mumlab;
+      gmumpdg[0] = reftrk.gmumpdg;
+      gmumlab[0] = reftrk.gmumlab;
+      refTree->GetEvent(theB->GetTrkNum(1));
+      pdg[1]    = reftrk.pdg;
+      mumpdg[1] = reftrk.mumpdg;
+      mumlab[1] = reftrk.mumlab;
+      gmumpdg[1] = reftrk.gmumpdg;
+      gmumlab[1] = reftrk.gmumlab;
+      delete refTree;
+    }
+    
+    theB->SetPdgCodes(pdg);
+    theB->SetMumPdgCodes(mumpdg);
+    theB->SetGMumPdgCodes(gmumpdg);
+
+    if(gmumpdg[0]==gmumpdg[1] &&              // Both GrandMothers are of the same sign
+       (TMath::Abs(gmumpdg[0])==521   || TMath::Abs(gmumpdg[0])==511   ||  // GrandMother Bplus/Bminus or B0/B0bar
+        TMath::Abs(gmumpdg[0])==523   || TMath::Abs(gmumpdg[0])==513   ||  // B0s/B0sbar
+        TMath::Abs(gmumpdg[0])==515   || TMath::Abs(gmumpdg[0])==525   ||  // 
+        TMath::Abs(gmumpdg[0])==531   || TMath::Abs(gmumpdg[0])==533   ||  // 
+        TMath::Abs(gmumpdg[0])==535   || TMath::Abs(gmumpdg[0])==541   ||  // 
+        TMath::Abs(gmumpdg[0])==543   || TMath::Abs(gmumpdg[0])==545   ||  // 
+        TMath::Abs(gmumpdg[0])==10521 || TMath::Abs(gmumpdg[0])==10511 ||  //   all possible 
+        TMath::Abs(gmumpdg[0])==10523 || TMath::Abs(gmumpdg[0])==10513 ||  //    B mesons
+        TMath::Abs(gmumpdg[0])==20523 || TMath::Abs(gmumpdg[0])==20513 ||  // 
+        TMath::Abs(gmumpdg[0])==10531 || TMath::Abs(gmumpdg[0])==10533 ||  // 
+        TMath::Abs(gmumpdg[0])==20533 || TMath::Abs(gmumpdg[0])==10541 ||  // 
+        TMath::Abs(gmumpdg[0])==20543 || TMath::Abs(gmumpdg[0])==10543 ||  // 
+        TMath::Abs(gmumpdg[0])==4122  || TMath::Abs(gmumpdg[0])==4222  ||  // All possible B baryons
+        TMath::Abs(gmumpdg[0])==4212  || TMath::Abs(gmumpdg[0])==4112  ||  // All possible B baryons
+        TMath::Abs(gmumpdg[0])==4224  || TMath::Abs(gmumpdg[0])==4214  ||  // All possible B baryons
+        TMath::Abs(gmumpdg[0])==4114  || TMath::Abs(gmumpdg[0])==4232  ||  // All possible B baryons
+        TMath::Abs(gmumpdg[0])==4132  || TMath::Abs(gmumpdg[0])==4322  ||  // All possible B baryons
+        TMath::Abs(gmumpdg[0])==4312  || TMath::Abs(gmumpdg[0])==4324  ||  // All possible B baryons
+        TMath::Abs(gmumpdg[0])==4314  || TMath::Abs(gmumpdg[0])==4332  ||  // All possible B baryons
+        TMath::Abs(gmumpdg[0])==4334  || TMath::Abs(gmumpdg[0])==4412  ||  // All possible B baryons
+        TMath::Abs(gmumpdg[0])==4422  || TMath::Abs(gmumpdg[0])==4414  ||  // All possible B baryons
+        TMath::Abs(gmumpdg[0])==4424  || TMath::Abs(gmumpdg[0])==4432  ||  // All possible B baryons
+        TMath::Abs(gmumpdg[0])==4434  || TMath::Abs(gmumpdg[0])==4444      // All possible B baryons
+       ) &&
+       mumpdg[0]==443 && mumpdg[1]== 443 && 
+       mumlab[0]==mumlab[1] && 
+       reftrk.mumprongs==2 &&
+       pdg[0]==-11 && pdg[1]==11  
+     ) theB->SetSignal();
+
+    else if (  // here consider the case of primary J/psi 
+      gmumlab[0]<0 && gmumlab[1]<0 && 
+       mumpdg[0]==443 && mumpdg[1]== 443 &&
+       mumlab[0]==mumlab[1] &&
+       reftrk.mumprongs==2 &&
+       pdg[0]==-11 && pdg[1]==11
+     ) theB->SetJpsiPrimary();
+
+    // if(!fOnlySignal || theB->IsSignal()) treeBout->Fill();  
+    
+ // write it out 1) always if you have not asked for onlySignal or OnlyPrimaryJpsi (or both)
+ //          or  2) if you have asked for Signal and it is Signal
+ //          or  3) if you have asked for Primary Jpsi and it is a Primary Jpsi
+    if ( (!fOnlySignal && !fOnlyPrimaryJpsi) || (fOnlySignal && theB->IsSignal()) 
+        || (fOnlyPrimaryJpsi && theB->IsJpsiPrimary()) ) treeBout->Fill();
+  }
+
+  delete refTree0;
+
+  refFile->Close();
+
+  return;
+}
+
+
+
+
+
+
diff --git a/PWG3/AliBtoJPSItoEleAnalysis.h b/PWG3/AliBtoJPSItoEleAnalysis.h
new file mode 100644 (file)
index 0000000..b534b15
--- /dev/null
@@ -0,0 +1,92 @@
+#ifndef AliBtoJPSItoEleAnalysis_H
+#define AliBtoJPSItoEleAnalysis_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                      Class AliBtoJPSItoEleAnalysis
+//             Reconstruction and analysis B -> J/\psi + X 
+//                                               |_> e+ e- 
+//      
+//         Origin: G.E Bruno    giuseppe.bruno@ba.infn.it                
+//  based on Class for charm golden channel (D0->Kpi)
+//-------------------------------------------------------------------------
+
+#include <TString.h>
+#include <TNamed.h>
+#include "AliESD.h"
+#include "AliRun.h"
+
+//-----------------------------------------------------------------------------
+class AliBtoJPSItoEleAnalysis : public TNamed {
+ public:
+  //
+  AliBtoJPSItoEleAnalysis();
+  virtual ~AliBtoJPSItoEleAnalysis();
+
+  void ApplySelection(const Char_t *inName="AliBtoJPSItoEle.root",
+                     const Char_t *outName="AliBtoJPSItoEle_sele.root") const;
+  void FindCandidates(Int_t evFirst=0,Int_t evLast=0,
+                     const Char_t *outName="AliBtoJPSItoEle.root");
+  void MakeTracksRefFile(AliRun *gAlice,Int_t evFirst=0,Int_t evLast=0) const;
+  void PrintStatus() const;
+  void SetVertexOnTheFly() { fVertexOnTheFly=kTRUE; }
+  void SetSimulation() { fSim=kTRUE; }
+  void SetOnlySignal() { fOnlySignal=kTRUE; }
+  void SetOnlyPrimaryJpsi() { fOnlyPrimaryJpsi=kTRUE; }
+  void SetOnlySignalAndPrimaryJpsi() { fOnlySignal=kTRUE; fOnlyPrimaryJpsi=kTRUE; }
+  void SetPtCut(Double_t pt=0.) { fPtCut=pt; }
+  void Setd0Cut(Double_t d0=0.) { fd0Cut=d0; } 
+  void SetPidCut(Double_t pid=0.) { fPidCut=pid; }
+  void SetMassCut(Double_t deltaM=1000.) { fMassCut=deltaM; }
+  void SetBCuts(Double_t cut0=1000.,Double_t cut1=100000.,
+                Double_t cut2=1.1,Double_t cut3=0.,Double_t cut4=0.,
+                Double_t cut5=100000.,Double_t cut6=100000.,
+                Double_t cut7=100000000.,Double_t cut8=-1.1); 
+  void SetBCuts(const Double_t cuts[9]); 
+  void SetPID(const Char_t * pid="TRDTPCparam") { fPID=pid; }
+  //
+ private:
+  //
+  Bool_t   fVertexOnTheFly; // flag for primary vertex reco on the fly
+  Bool_t   fSim;            // flag for the analysis of simulated events
+  Bool_t   fOnlySignal;     // write to file only signal candidates (for sim)
+  Bool_t   fOnlyPrimaryJpsi;// write to file only primary Jpsi candidates (for sim)
+  TString  fPID;           // PID scheme
+
+  Double_t fV1[3]; // primary vertex position (in cm)
+  Double_t fPtCut;   // minimum track pt (in GeV/c)
+  Double_t fd0Cut;   // minimum track |rphi impact parameter| (in micron) 
+  Double_t fMassCut; // maximum of |InvMass-M(J/Psi)| (in GeV)
+  Double_t fPidCut; // min. pid probability as an electron
+  Double_t fBCuts[9]; // cuts on b candidates (see SetBCuts())
+                       // (to be passed to function AliBtoJPSItoEle::Select())
+                       // 0 = inv. mass half width [GeV]   
+                       // 1 = dca [micron]
+                       // 2 = cosThetaStar 
+                       // 3 = pTP [GeV/c] (positron)
+                       // 4 = pTN [GeV/c] (electron)
+                       // 5 = d0P [micron]   upper limit!
+                       // 6 = d0N [micron]  upper limit!
+                       // 7 = d0d0 [micron^2]
+                       // 8 = cosThetaPoint
+
+  //
+  Double_t CalculateTOFmass(Double_t mom,Double_t length,Double_t time) const;
+  Bool_t   SelectInvMass(const Double_t p[6]) const;
+  void     SelectTracks(AliESD *event,
+                       TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
+                       TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const;
+  void     SetVertex1(Double_t x=0.,Double_t y=0.,Double_t z=0.) 
+    { fV1[0]=x;fV1[1]=y;fV1[2]=z; }
+  void     SimulationInfo(TTree *treeBin,TTree *treeBout) const;
+  Bool_t   SingleTrkCuts(const AliESDtrack& trk, Double_t b) const;
+  //
+  ClassDef(AliBtoJPSItoEleAnalysis,1)  // Reconstruction of B->JPSI-> e+e- candidates class
+};
+
+
+#endif
+
+
+
diff --git a/PWG3/AliBtoJPSItoEleReco_all.C b/PWG3/AliBtoJPSItoEleReco_all.C
new file mode 100644 (file)
index 0000000..e508fe0
--- /dev/null
@@ -0,0 +1,94 @@
+//--------------------------------------------------------------------------
+// Test macro for reconstruction and analysis of B->J/psi*X
+//                                                  |_>e+e-
+//
+//     Giuseppe Bruno, giuseppe.bruno@ba.infn.it
+//      based on the for charm golden channel (D0->Kpi)
+//--------------------------------------------------------------------------
+
+void AliBtoJPSItoEleReco_all() {
+  
+  gSystem->Load("libAOD.so");
+  gSystem->Load("libPWG3base.so");
+
+  //==============  R E C O N S T R U C T I O N ==============================
+
+  Int_t evFirst = 0;
+  Int_t evLast  = 99;
+
+  // Get field from galice.root
+  if (gAlice) {
+    delete gAlice->GetRunLoader();
+    delete gAlice; 
+    gAlice=0;
+  }  
+  AliRunLoader *rl = AliRunLoader::Open("galice.root");
+  if (rl == 0x0) {
+    cerr<<"Can not open session"<<endl;
+    return;
+  }
+  Int_t retval = rl->LoadgAlice();
+  if (retval) {
+    cerr<<"LoadgAlice returned error"<<endl;
+    delete rl;
+    return;
+  }
+  gAlice=rl->GetAliRun();
+  AliMagF *fiel = (AliMagF*)gAlice->Field();
+  // Set the conversion constant between curvature and Pt
+  AliTracker::SetFieldMap(fiel,kTRUE);
+
+  AliBtoJPSItoEleAnalysis *analysis = new AliBtoJPSItoEleAnalysis();
+  // set simulation to take info on PDG codes from kine tree
+  analysis->SetSimulation();
+  rl->LoadKinematics(); 
+  analysis->MakeTracksRefFile(gAlice,evFirst,evLast);
+  // cout << "ci arrivo 2" << endl;
+  //--- set this is you want only signal candidates in output file
+   analysis->SetOnlySignal();
+  //--- set this if you want to compute primary vertex candidate by candidate using 
+  //    other tracks in the event (for pp, broad interaction region);
+  //    it is time-consuming procedure, so it can be done after a 
+  //    preselection on invariant mass
+  //analysis->SetVertexOnTheFly();
+  //analysis->SetMassCut(0.1); // GeV
+  //--- set single-track preselections
+  analysis->SetPtCut(0.); // GeV
+  analysis->Setd0Cut(0.); // micron
+  //--- set cuts on candidates to be written to file
+  //    (see AliBtoJPSItoEleAnalysis.h for a description and for the defaults)
+  //analysis->SetBCuts(0.1,1000.,1.1,0.,0.,10000.,10000.,0.,.5);
+  analysis->SetBCuts();
+
+  //--- check the settings
+  analysis->PrintStatus();
+
+  analysis->FindCandidates(evFirst,evLast,"AliBtoJPSItoEle_all.root");
+  delete analysis;
+
+  return;
+}
+//==========================================================================
+void AliBtoJPSItoEleSele() {  
+
+  gSystem->Load("libAOD.so");
+  gSystem->Load("libPWG3base.so");
+
+  //========================  S E L E C T I O N ============================
+
+  AliBtoJPSItoEleAnalysis *analysis = new AliBtoJPSItoEleAnalysis();
+  analysis->SetSimulation();
+  analysis->SetOnlySignal();
+  analysis->SetBCuts(0.1,1000.,1.1,0.,0.,10000.,10000.,0.,.5);
+  analysis->ApplySelection("AliBtoJPSItoEle.root","AliBtoJPSItoEle_sele.root");
+  delete analysis;
+
+  return;
+}
+//==========================================================================
+
+
+
+
+
diff --git a/PWG3/AliPrimJPSItoEleReco_all.C b/PWG3/AliPrimJPSItoEleReco_all.C
new file mode 100644 (file)
index 0000000..4ae8300
--- /dev/null
@@ -0,0 +1,96 @@
+//--------------------------------------------------------------------------
+// Test macro for reconstruction and analysis of Primary J/psi into e+e-
+//
+//     Giuseppe Bruno, giuseppe.bruno@ba.infn.it
+//      based on the for charm golden channel (D0->Kpi)
+//--------------------------------------------------------------------------
+
+void AliPrimJPSItoEleReco_all() {
+  
+  gSystem->Load("libAOD.so");
+  gSystem->Load("libPWG3base.so");
+
+  //==============  R E C O N S T R U C T I O N ==============================
+
+  Int_t evFirst = 0;
+  Int_t evLast  = 99;
+
+  // Get field from galice.root
+  if (gAlice) {
+    delete gAlice->GetRunLoader();
+    delete gAlice; 
+    gAlice=0;
+  }  
+  AliRunLoader *rl = AliRunLoader::Open("galice.root");
+  if (rl == 0x0) {
+    cerr<<"Can not open session"<<endl;
+    return;
+  }
+  Int_t retval = rl->LoadgAlice();
+  if (retval) {
+    cerr<<"LoadgAlice returned error"<<endl;
+    delete rl;
+    return;
+  }
+  gAlice=rl->GetAliRun();
+  AliMagF *fiel = (AliMagF*)gAlice->Field();
+  // Set the conversion constant between curvature and Pt
+  AliTracker::SetFieldMap(fiel,kTRUE);
+
+  AliBtoJPSItoEleAnalysis *analysis = new AliBtoJPSItoEleAnalysis();
+  // set simulation to take info on PDG codes from kine tree
+  analysis->SetSimulation();
+  rl->LoadKinematics(); 
+  analysis->MakeTracksRefFile(gAlice,evFirst,evLast);
+  // cout << "ci arrivo 2" << endl;
+  //--- set this is you want only signal candidates in output file
+   //analysis->SetOnlySignal();
+   analysis->SetOnlyPrimaryJpsi();
+   //analysis->SetOnlySignalAndPrimaryJpsi();
+  //--- set this if you want to compute primary vertex candidate by candidate using 
+  //    other tracks in the event (for pp, broad interaction region);
+  //    it is time-consuming procedure, so it can be done after a 
+  //    preselection on invariant mass
+  //analysis->SetVertexOnTheFly();
+  //analysis->SetMassCut(0.1); // GeV
+  //--- set single-track preselections
+  analysis->SetPtCut(0.); // GeV
+  analysis->Setd0Cut(0.); // micron
+  //--- set cuts on candidates to be written to file
+  //    (see AliBtoJPSItoEleAnalysis.h for a description and for the defaults)
+  //analysis->SetBCuts(0.1,1000.,1.1,0.,0.,10000.,10000.,0.,.5);
+  analysis->SetBCuts();
+
+  //--- check the settings
+  analysis->PrintStatus();
+
+  analysis->FindCandidates(evFirst,evLast,"AliJPSItoEle_all.root");
+  delete analysis;
+
+  return;
+}
+//==========================================================================
+void AliBtoJPSItoEleSele() {  
+
+  gSystem->Load("libAOD.so");
+  gSystem->Load("libPWG3base.so");
+
+  //========================  S E L E C T I O N ============================
+
+  AliBtoJPSItoEleAnalysis *analysis = new AliBtoJPSItoEleAnalysis();
+  analysis->SetSimulation();
+  //analysis->SetOnlySignal();
+  analysis->SetOnlyPrimaryJpsi();
+  analysis->SetBCuts(0.1,1000.,1.1,0.,0.,10000.,10000.,0.,.5);
+  analysis->ApplySelection("AliBtoJPSItoEle.root","AliBtoJPSItoEle_sele.root");
+  delete analysis;
+
+  return;
+}
+//==========================================================================
+
+
+
+
+