]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
1) possibility to use KF vertexer(Giuseppe)
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Nov 2007 13:21:31 +0000 (13:21 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Nov 2007 13:21:31 +0000 (13:21 +0000)
PWG3/AliBtoJPSItoEle.cxx
PWG3/AliBtoJPSItoEle.h
PWG3/AliBtoJPSItoEleAnalysis.cxx
PWG3/AliBtoJPSItoEleAnalysis.h

index 842e154bcfb6e28aea35895514272a6ce7d0f23d..a17583a2fa5f026d1d97411826e1f959c2cb2172 100644 (file)
 #include <TCanvas.h>
 #include <TPaveLabel.h>
 #include <TVector3.h>
+#include <TString.h>
 
 #include "AliBtoJPSItoEle.h"
 
 ClassImp(AliBtoJPSItoEle)
 
 //----------------------------------------------------------------------------
-AliBtoJPSItoEle::AliBtoJPSItoEle():
-fSignal(kFALSE),
-fJpsiPrimary(kFALSE),
-fEvent(0),
-fV1x(0.),
-fV1y(0.),
-fV1z(0.),
-fV2x(0.),
-fV2y(0.),
-fV2z(0.),
-fDCA(0.),
-fWgtJPsi(0.)
-{
+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.;
@@ -77,29 +79,31 @@ fWgtJPsi(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]):
-fSignal(kFALSE),
-fJpsiPrimary(kFALSE),
-fEvent(ev),
-fV1x(v1[0]),
-fV1y(v1[1]),
-fV1z(v1[2]),
-fV2x(v2[0]),
-fV2y(v2[1]),
-fV2z(v2[2]),
-fDCA(dca),
-fWgtJPsi(0.)
-{
+                      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];
@@ -123,9 +127,15 @@ fWgtJPsi(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
@@ -380,6 +390,14 @@ void AliBtoJPSItoEle::SetPIDresponse(Double_t resp0[5],Double_t resp1[5]) {
   return;
 } 
 //-----------------------------------------------------------------------------
+Double_t AliBtoJPSItoEle::GetTagEl(Int_t child) const {
+  // Get tag probabilities for electrons
+                                                                                                                          
+  if(child!=0 && child !=1) return -1;
+
+  return fTagEl[child];
+}
+//-----------------------------------------------------------------------------
 void AliBtoJPSItoEle::GetPIDresponse(Double_t resp0[5],Double_t resp1[5]) const {
   // Get combined PID detector response probabilities
                                                                                                                           
@@ -412,10 +430,13 @@ Double_t AliBtoJPSItoEle::TRDTPCCombinedPIDParametrization(Double_t p) const {
   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.;
+// Better estimation based on TRD test beam (as presented by Andrea at Munster)
+//  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.;
+// From Silvia MASCIOCCHI (Darmstadt) at PWG3 AliceWeek October 2007
+  if (p<10.) value*=(0.44-0.06*p+0.031*p*p-0.0008*p*p*p)/100.;
+  if (p>10.) value*=(-0.67+0.28*p)/100.;
 
   return value;
 }
@@ -424,4 +445,3 @@ Double_t AliBtoJPSItoEle::TRDTPCCombinedPIDParametrization(Double_t p) const {
 
 
 
-
index fb6dcb9481df2ef4ad943e5cadcda9a16b7615a3..5272a62c819f5719ebc34b87d45819d144c6ed8a 100644 (file)
@@ -36,6 +36,7 @@ class AliBtoJPSItoEle : public TObject {
             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");
@@ -57,6 +58,7 @@ class AliBtoJPSItoEle : public TObject {
   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; 
+  Double_t GetTagEl(Int_t child) const; 
   void     GetWgts(Double_t &WgtJPsi) const;
   void     GetPrimaryVtx(Double_t vtx[3]) const 
     { vtx[0]=fV1x; vtx[1]=fV1y; vtx[2]=fV1z; }
@@ -147,4 +149,3 @@ class AliBtoJPSItoEle : public TObject {
 
 
 
-
index 81cf2e4917286bad68c94eabe34b64ed94136496..70771f07c3012579def100139395b98ecab3e85f 100644 (file)
@@ -37,6 +37,9 @@
 #include "AliBtoJPSItoEle.h"
 #include "AliBtoJPSItoEleAnalysis.h"
 #include "AliLog.h"
+#include "AliKFParticleBase.h"
+#include "AliKFParticle.h"
+#include "AliKFVertex.h"
 
 typedef struct {
   Int_t lab;
@@ -62,7 +65,10 @@ fPID("TRDTPCparam"),
 fPtCut(0.),
 fd0Cut(0.), 
 fMassCut(1000.),
-fPidCut(0.)
+fPidCut(0.),
+//fKFPrimVertex(kFALSE),
+//fKFTopConstr(kFALSE),
+fKFSecondVertex(kFALSE)
 {
   // Default constructor
 
@@ -195,6 +201,27 @@ void AliBtoJPSItoEleAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
   }
   event->ReadFromTree(tree);
 
+/*  if (fKFPrimVertex)
+  AliRunLoader* runLoader = 0;
+  {
+    if (gAlice) {
+      delete gAlice->GetRunLoader();
+      delete gAlice;
+      gAlice=0;
+    }
+    runLoader = AliRunLoader::Open(galName.Data());
+    if (runLoader == 0x0) {
+      cerr<<"Can not open session"<<endl;
+      return;
+    }
+    cout << "Ok open galice.root" << endl;
+    runLoader->LoadgAlice();
+
+    gAlice = runLoader->GetAliRun();
+    runLoader->LoadKinematics();
+    runLoader->LoadHeader();
+  } */
+
   // loop on events in file
   for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
     if(iEvent > evLast) break;
@@ -230,6 +257,87 @@ void AliBtoJPSItoEleAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
 
     nBrec1ev = 0;
 
+/*
+//===================== PRIMARY VERTEX USING KF METHODS ==========================//
+
+if (fKFPrimVertex)
+{
+  AliStack* stack = runLoader->Stack();
+
+  class TESDTrackInfo
+    {
+    public:
+    TESDTrackInfo(){}
+    AliKFParticle fParticle; // KFParticle constructed from ESD track
+    //Bool_t fPrimUsedFlag;    // flag says that the particle was used for primary vertex fit
+    Bool_t fOK;              // is the track good enough
+    Int_t mcPDG;             // Monte Carlo PDG code of the particle
+    Int_t mcMotherID;        // Monte Carlo ID of its mother
+    };
+
+ // TESDTrackInfo ESDTrackInfo[trkEntries];
+  TESDTrackInfo ESDTrackInfo[1000];
+  for (Int_t iTr=0; iTr<trkEntries; iTr++)
+    {
+      TESDTrackInfo &info = ESDTrackInfo[iTr];
+      info.fOK = 0;
+      //info.fPrimUsedFlag = 0;
+      info.mcPDG = -1;
+      info.mcMotherID = -1;
+
+      // track quality check
+
+      AliESDtrack *pTrack = event->GetTrack(iTr);
+      if( !pTrack  ) continue;
+      if (pTrack->GetKinkIndex(0)>0) continue;
+      if ( !( pTrack->GetStatus()&AliESDtrack::kITSrefit ) ) continue;
+      //Int_t indi[12];
+      //if( pTrack->GetITSclusters(indi) <5 ) continue;
+      //Int_t PDG = ( pTrack->GetSigned1Pt() <0 ) ?321 :211;
+
+      // take MC PDG
+
+      Int_t mcID = TMath::Abs(pTrack->GetLabel());
+      TParticle * part = stack->Particle(TMath::Abs(mcID));
+      info.mcPDG = part->GetPdgCode();
+      Int_t PDG = info.mcPDG;
+      if( mcID>=0 ) info.mcMotherID = part->GetFirstMother();
+
+
+      // Construct KFParticle for the track
+
+      info.fParticle = AliKFParticle( *pTrack, PDG );
+      info.fOK = 1;
+    }
+
+    // Find event primary vertex with KF methods
+
+   AliKFVertex primVtx;
+    {
+   // const AliKFParticle * vSelected[trkEntries]; // Selected particles for vertex fit
+   // Int_t vIndex[trkEntries];                    // Indices of selected particles
+   // Bool_t vFlag[trkEntries];                    // Flags returned by the vertex finder
+   const AliKFParticle * vSelected[1000]; // Selected particles for vertex fit
+   Int_t vIndex[1000];                    // Indices of selected particles
+   Bool_t vFlag[1000];                    // Flags returned by the vertex finder
+
+    Int_t nSelected = 0;
+    for( Int_t i = 0; i<trkEntries; i++){
+      if(ESDTrackInfo[i].fOK ){
+        vSelected[nSelected] = &(ESDTrackInfo[i].fParticle);
+        vIndex[nSelected] = i;
+        nSelected++;
+      }
+    }
+    primVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, 3. );
+    //for( Int_t i = 0; i<nSelected; i++){
+    //  if( vFlag[i] ) ESDTrackInfo[vIndex[i]].fPrimUsedFlag = 1;
+    //}
+    if( primVtx.GetNDF() <1 ) return; // Less then two tracks in primary vertex
+    }
+
+}*/
+
     // 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);
@@ -257,12 +365,25 @@ void AliBtoJPSItoEleAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
        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);
+
+        if (fKFSecondVertex){
+          //Define the AliKFParticle Objects
+          AliKFParticle trackP = AliKFParticle(pt,-11);
+          AliKFParticle trackN = AliKFParticle(nt,11);
+          //Construct the V0like mother
+          AliKFParticle V0(trackP,trackN);
+          //Get global position of the secondary vertex using KF methods
+          v2[0] = V0.GetX();
+          v2[1] = V0.GetY();
+          v2[2] = V0.GetZ();
+          mom[0] = trackP.GetPx(); mom[1] = trackP.GetPy(); mom[2] = trackP.GetPz();
+          mom[3] = trackN.GetPx(); mom[4] = trackN.GetPy(); mom[5] = trackN.GetPz();
+        }else{
+          //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]);
+        }
        goodVtx1 = kTRUE;
        
        // no vertexing if DeltaMass > fMassCut 
@@ -283,6 +404,28 @@ void AliBtoJPSItoEleAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
          }
        }
        
+/*      if (fKFSecondVertex&&fKFTopConstr&&fKFPrimVertex){
+//====================== TOPOLOGICAL CONSTRAINT !!=====================================
+        //Primary vertex constructed from ESD using KF methods!!!
+        AliKFVertex primVtxCopy(*(event->GetPrimaryVertex()));
+
+        //Subtract Daughters from primary vertex
+        primVtxCopy -= trackP;
+        primVtxCopy -= trackN;
+
+        //Add V0 to the vertex in order to improve primary vertex resolution
+        primVtxCopy += V0;
+
+        //Set production vertex for V0
+        V0.SetProductionVertex(primVtxCopy);
+
+        //Recalculate primary vertex
+        fV1[0] = primVtxCopy.GetX();
+        fV1[1] = primVtxCopy.GetY();
+        fV1[2] = primVtxCopy.GetZ();
+//=====================================================================================
+        }*/
+
        // 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);
@@ -693,11 +836,18 @@ void AliBtoJPSItoEleAnalysis::SimulationInfo(TTree *treeBin,TTree *treeBout) con
      ) theB->SetSignal();
 
     else if (  // here consider the case of primary J/psi 
-      gmumlab[0]<0 && gmumlab[1]<0 && 
        mumpdg[0]==443 && mumpdg[1]== 443 &&
+       pdg[0]==-11 && pdg[1]==11 &&
        mumlab[0]==mumlab[1] &&
        reftrk.mumprongs==2 &&
-       pdg[0]==-11 && pdg[1]==11
+       ( gmumlab[0]<0                     ||   // really primary J/psi (without family. e.g. from Cocktail)
+         TMath::Abs(gmumpdg[0])==100443   ||   // from Psi(2S)
+         TMath::Abs(gmumpdg[0])==10441    ||   // from Csi_c0(1P)
+         TMath::Abs(gmumpdg[0])==20443    ||   // from Csi_c1(1P)
+         TMath::Abs(gmumpdg[0])==10443    ||   // from h_c(1P)
+         TMath::Abs(gmumpdg[0])==445      ||   // from Csi_c2(1P)
+         (gmumpdg[0]>=81 && gmumpdg[0]<=100)   // this is for MC internal use (e.g. J/psi from string)
+       )
      ) theB->SetJpsiPrimary();
 
     // if(!fOnlySignal || theB->IsSignal()) treeBout->Fill();  
@@ -720,4 +870,3 @@ void AliBtoJPSItoEleAnalysis::SimulationInfo(TTree *treeBin,TTree *treeBout) con
 
 
 
-
index 5c5c87104e536fa0e1a4512468ca7bc8fbff2151..c550745fa9b941e903739a110e73dda4b9307b35 100644 (file)
@@ -45,6 +45,10 @@ class AliBtoJPSItoEleAnalysis : public TNamed {
                 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; }
+//  void SetKFPrimVertex() { fKFPrimVertex=kTRUE; }          //new setter
+   void SetKFSecondVertex() {fKFSecondVertex=kTRUE;}         //new setter
+   void UnSetKFSecondVertex() {fKFSecondVertex=kFALSE;}      //new setter
+//  void SetKFTopConstr() { fKFTopConstr=kTRUE; }            //new setter
   //
  private:
   //
@@ -59,6 +63,9 @@ class AliBtoJPSItoEleAnalysis : public TNamed {
   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
+  Bool_t   fKFSecondVertex;  // flag for Kalmann Filter reco of secondary vertex
+//  Bool_t   fKFTopConstr;     // flag for Kalmann Filter topological constraint in primary vtx reco
+//  Bool_t   fKFPrimVertex;    // flag for Kalmann Filter reco of primary vertex
   Double_t fBCuts[9]; // cuts on b candidates (see SetBCuts())
                        // (to be passed to function AliBtoJPSItoEle::Select())
                        // 0 = inv. mass half width [GeV]   
@@ -82,11 +89,10 @@ class AliBtoJPSItoEleAnalysis : public TNamed {
   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
+  ClassDef(AliBtoJPSItoEleAnalysis,2)  // Reconstruction of B->JPSI-> e+e- candidates class
 };
 
 
 #endif
 
 
-