added more options to run D0 finder from offline or online data (Gaute)
authorrichterm <richterm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 26 Nov 2009 08:24:17 +0000 (08:24 +0000)
committerrichterm <richterm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 26 Nov 2009 08:24:17 +0000 (08:24 +0000)
HLT/trigger/AliHLTD0Trigger.cxx
HLT/trigger/AliHLTD0Trigger.h
HLT/trigger/AliHLTD0toKpi.cxx
HLT/trigger/AliHLTD0toKpi.h

index dee05d6..97b3363 100644 (file)
@@ -28,7 +28,6 @@
 
 #include "AliHLTD0Trigger.h"
 #include "AliESDEvent.h"
-#include "AliESDtrack.h"
 #include "AliESDv0.h"
 #include "AliHLTTriggerDecision.h"
 #include "AliHLTDomainEntry.h"
@@ -39,6 +38,7 @@
 #include "AliESDVertex.h"
 #include "TH1F.h"
 #include "AliHLTD0toKpi.h"
+#include "AliAODVertex.h"
 
 /** ROOT macro for the implementation of ROOT specific class methods */
 ClassImp(AliHLTD0Trigger)
@@ -52,9 +52,17 @@ AliHLTD0Trigger::AliHLTD0Trigger()
   , fd0(0.0) 
   , fd0d0(0.0)
   , fcosPoint(0.0)
-  ,fplothisto(false)
+  , fplothisto(false)
+  , fUseV0(false)
+  , mD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
   , fD0mass(NULL)
+  , fPos()
+  , fNeg()
+  , fd0calc(NULL)                  
+  , ftwoTrackArray(NULL)
+  , fTotalD0(0)
 {
+  
   // see header file for class documentation
   // or
   // refer to README to build package
@@ -67,10 +75,12 @@ const char* AliHLTD0Trigger::fgkOCDBEntry="HLT/ConfigHLT/D0Trigger";
 
 AliHLTD0Trigger::~AliHLTD0Trigger()
 {
+  //if(fd0calc){delete fd0calc;}  
+  //if(fD0mass){delete fD0mass;}
+  //if(ftwoTrackArray){delete ftwoTrackArray;}
   // see header file for class documentation
-  //delete fD0mass;
 }
-
+  
 const char* AliHLTD0Trigger::GetTriggerName() const
 {
   // see header file for class documentation
@@ -85,73 +95,42 @@ AliHLTComponent* AliHLTD0Trigger::Spawn()
 
 int AliHLTD0Trigger::DoTrigger()
 {
-  // see header file for class documentation
-  int iResult=0;
-  int nD0=0;
-  Double_t D0,D0bar; 
-  Double_t d0[2];
-  Double_t svpos[3];
-  Double_t pvpos[3];
-  TString description;
-  AliHLTD0toKpi *d0calc = new AliHLTD0toKpi();
-  
-  for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
+  // -- Iterator over Data Blocks --
+  //const AliHLTComponentBlockData* iter = NULL;
 
-    AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
-    event->GetStdContent();
-    Int_t nV0 = event->GetNumberOfV0s();
-  
-    Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
-
-    for (Int_t iv=0; iv<nV0; iv++) {
+  if (!IsDataEvent()) return 0;
 
-      AliESDtrack *tN=event->GetTrack( event->GetV0(iv)->GetNindex());
-      AliESDtrack *tP=event->GetTrack( event->GetV0(iv)->GetPindex());      
-      
-      if(tN->Pt()<fPtMin && tP->Pt()<fPtMin){continue;}
-      
-      const AliESDVertex* pv = event->GetPrimaryVertexTracks();
-      Double_t field = event->GetMagneticField();
-
-      pv->GetXYZ(pvpos);
-
-      d0[0] =  10000.*tP->GetD(pvpos[0],pvpos[1],field);
-      d0[1] = -10000.*tN->GetD(pvpos[0],pvpos[1],field);
-
-      if(d0[0]<fd0 && d0[0]<fd0){continue;} // make sure < or>
-
-      event->GetV0(iv)->GetXYZ(svpos[0],svpos[1],svpos[2]);
+  Int_t nD0=0;
+  TString description;
 
-      if(!tN->PropagateTo(svpos[0],field) && !tP->PropagateTo(svpos[0],field)){
-       HLTInfo("Tracks could not be propagated to secondary vertex");
-       continue;
-      }
-      
-      Double_t tmp1, tmp2;
-      if(tN->GetDCA(tP,field,tmp1,tmp2) > fdca){continue;}
+    HLTDebug("Cuts: -pt:%f -dca:%f -invmass:%f -costhetastar:%f -d0:%f -d0d0:%f -cospoint:%f",fPtMin,fdca,finvMass,fcosThetaStar,fd0,fd0d0,fcosPoint);
 
-      if((d0calc->InvMass(tN,tP) - mD0PDG) > finvMass && (d0calc->InvMass(tP,tN) - mD0PDG) > finvMass){continue;}
-      d0calc->cosThetaStar(tN,tP,D0,D0bar);
-      if(D0 > fcosThetaStar && D0bar > fcosThetaStar){continue;}
-      if((d0[0]*d0[1]) > fd0d0){continue;}
-      if(d0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
-      
-      nD0++;
-      if((d0calc->InvMass(tN,tP) - mD0PDG) > finvMass){
-       fD0mass->Fill(d0calc->InvMass(tN,tP));
-      }
-      else{
-       fD0mass->Fill(d0calc->InvMass(tP,tN));
-      }
+  for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {   
+    if(fUseV0){
+      nD0=RecV0(iter);
+    }
+    else{
+      nD0=RecESDTracks(iter);
     }
-
   }
-  HLTWarning("Number of D0 found: %d",nD0);
+
+  for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeTrack); iter != NULL; iter = GetNextInputObject() ) { 
+    nD0=RecBarrelTracks(iter);
+  }    
+
+  fTotalD0+=nD0;
+   
+  ftwoTrackArray->Clear();
+  fPos.clear();
+  fNeg.clear();
+  
+  HLTDebug("Number of D0 found: %d",nD0);
+  HLTDebug("Total Number of D0 found: %d",fTotalD0);
 
   if(fplothisto){PushBack( (TObject*) fD0mass, kAliHLTDataTypeHistogram,0);}
   
-  if (iResult>=0) {
+  //if (iResult>=0) {
+  if (1) {
    
     if (nD0>=1) {
       description.Form("Event contains %d D0(s)", nD0);
@@ -178,14 +157,19 @@ int AliHLTD0Trigger::DoTrigger()
   }
   SetDescription(description.Data());
   TriggerEvent(false);
-  return iResult;
+  //return iResult;
+  return 0;
 }
 
 int AliHLTD0Trigger::DoInit(int argc, const char** argv)
 {
+  
+  fd0calc = new AliHLTD0toKpi();
+  ftwoTrackArray = new TObjArray(2);
+
   fplothisto=false;
   // see header file for class documentation
-  fD0mass = new TH1F("hD0mass","HLT:  D0 inv mass",500,0,50);
+  fD0mass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
   // first configure the default
   int iResult=0;
   if (iResult>=0) iResult=ConfigureFromCDBTObjString(fgkOCDBEntry);
@@ -199,7 +183,9 @@ int AliHLTD0Trigger::DoInit(int argc, const char** argv)
 int AliHLTD0Trigger::DoDeinit()
 {
   // see header file for class documentation
-  delete fD0mass;
+  if(fd0calc){delete fd0calc;}  
+  if(fD0mass){delete fD0mass;}
+  if(ftwoTrackArray){delete ftwoTrackArray;}
   return 0;
 }
 
@@ -274,11 +260,153 @@ int AliHLTD0Trigger::ScanConfigurationArgument(int argc, const char** argv)
     return 2;
   }
   if (argument.CompareTo("-plothistogram")==0) {
-    if (++i>=argc) return -EPROTO;
-    argument=argv[i];
     fplothisto=true;
-    return 2;
+    return 1;
+  }
+  if (argument.CompareTo("-usev0")==0) {
+    fUseV0=true;
+    return 1;
   }
+
   // unknown argument
   return -EINVAL;
 }
+
+void AliHLTD0Trigger::SingleTrackSelect(AliESDtrack* t, Double_t b,Double_t* pv){
+  // Offline har || på disse kuttene på de to henfallsproduktene 
+    
+  if(t->Pt()<fPtMin){return;}
+  if(TMath::Abs(t->GetD(pv[0],pv[1],b)) > fd0){return;}
+
+  if(t->Charge()>0){
+    fPos.push_back(t);
+  }
+  else{
+    fNeg.push_back(t);
+  }
+}
+
+Int_t AliHLTD0Trigger::RecESDTracks(const TObject* iter){
+  
+  int nD0=0;
+  Double_t D0,D0bar,xdummy,ydummy; 
+  Double_t d0[2];
+  Double_t svpos[3];
+  Double_t pvpos[3];
+  
+  AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
+  event->GetStdContent();
+  Double_t field = event->GetMagneticField();
+  const AliESDVertex* pv = event->GetPrimaryVertexTracks();
+  pv->GetXYZ(pvpos);
+  
+  for(Int_t it=0;it<event->GetNumberOfTracks();it++){
+    SingleTrackSelect(event->GetTrack(it),field,pvpos);
+  }
+  
+  for(Int_t ip=0;ip<fPos.size();ip++){
+    AliESDtrack *tP=fPos[ip];
+    for(Int_t in=0;in<fNeg.size();in++){
+      AliESDtrack *tN=fNeg[in];
+      
+      tP->PropagateToDCA(pv,field,kVeryBig);  //do I need this??????
+      tN->PropagateToDCA(pv,field,kVeryBig);
+      
+      Double_t dcatPtN = tP->GetDCA(tN,field,xdummy,ydummy);
+      if(dcatPtN>fdca) { continue; }
+      
+      ftwoTrackArray->AddAt(fPos[ip],0);
+      ftwoTrackArray->AddAt(fNeg[in],1);
+      AliAODVertex *vertexp1n1 = fd0calc->ReconstructSecondaryVertex(ftwoTrackArray,field,pv);
+      if(!vertexp1n1) { 
+       ftwoTrackArray->Clear();
+       continue; 
+      }
+      
+      vertexp1n1->GetXYZ(svpos);
+      
+      tP->PropagateToDCA(vertexp1n1,field,kVeryBig); 
+      tN->PropagateToDCA(vertexp1n1,field,kVeryBig);
+      
+      if((TMath::Abs(fd0calc->InvMass(tN,tP)-mD0PDG)) > finvMass && TMath::Abs((fd0calc->InvMass(tP,tN))-mD0PDG) > finvMass){continue;}
+      fd0calc->cosThetaStar(tN,tP,D0,D0bar);
+      if(TMath::Abs(D0) > fcosThetaStar && TMath::Abs(D0bar) > fcosThetaStar){continue;}
+      d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
+      d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
+      if((d0[0]*d0[1]) > fd0d0){continue;}
+      if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
+      
+      if(fplothisto){
+       if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass){
+         fD0mass->Fill(fd0calc->InvMass(tN,tP));
+       }
+       else{
+         fD0mass->Fill(fd0calc->InvMass(tP,tN));
+       }
+      }
+      nD0++;
+      delete vertexp1n1;
+    }
+  }
+  return nD0;
+}
+
+Int_t AliHLTD0Trigger::RecV0(const TObject* iter){
+  int nD0=0;
+  Double_t d0[2];
+  Double_t D0,D0bar; 
+  Double_t svpos[3];
+  Double_t pvpos[3];
+  
+  AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
+  event->GetStdContent();
+  Int_t nV0 = event->GetNumberOfV0s();
+  Double_t field = event->GetMagneticField();
+  const AliESDVertex* pv = event->GetPrimaryVertexTracks();
+  pv->GetXYZ(pvpos);
+
+  for (Int_t iv=0; iv<nV0; iv++) {
+      
+      AliESDtrack *tN=event->GetTrack( event->GetV0(iv)->GetNindex());
+      AliESDtrack *tP=event->GetTrack( event->GetV0(iv)->GetPindex());      
+      
+      if(tN->Pt()<fPtMin && tP->Pt()<fPtMin){continue;} //||????
+      
+      d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
+      d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
+      
+      if(d0[0]>fd0 && d0[0]>fd0){continue;} // make sure < or>
+      
+      event->GetV0(iv)->GetXYZ(svpos[0],svpos[1],svpos[2]);
+      
+      if(!tN->PropagateTo(svpos[0],field) && !tP->PropagateTo(svpos[0],field)){
+       HLTInfo("Tracks could not be propagated to secondary vertex");
+       continue;
+      }
+      
+      Double_t tmp1, tmp2;
+      if(tN->GetDCA(tP,field,tmp1,tmp2) > fdca){continue;}
+      
+      if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass && (fd0calc->InvMass(tP,tN) - mD0PDG) > finvMass){continue;}
+      fd0calc->cosThetaStar(tN,tP,D0,D0bar);
+      if(D0 > fcosThetaStar && D0bar > fcosThetaStar){continue;}
+      if((d0[0]*d0[1]) > fd0d0){continue;}
+      if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
+      
+      nD0++;
+      if((fd0calc->InvMass(tN,tP) - mD0PDG) > finvMass){
+       fD0mass->Fill(fd0calc->InvMass(tN,tP));
+      }
+      else{
+       fD0mass->Fill(fd0calc->InvMass(tP,tN));
+      }
+  }
+  return nD0;
+  
+}
+
+Int_t AliHLTD0Trigger::RecBarrelTracks(const TObject* iter){
+
+  return 0;
+
+}
index baf5bdb..f431733 100644 (file)
 /// @brief  HLT trigger component for D0->Kpi
 
 #include "AliHLTTrigger.h"
+#include <vector>
+#include "AliESDtrack.h"
+
 class TH1F;
+class TObjArray;
+class AliHLTD0toKpi;
 
 /**
  * @class  AliHLTD0Trigger
@@ -95,6 +100,11 @@ class AliHLTD0Trigger : public AliHLTTrigger
   /// inherited from AliHLTTrigger: calculate the trigger
   virtual int DoTrigger();
 
+  void SingleTrackSelect(AliESDtrack*,Double_t,Double_t*);
+  Int_t RecESDTracks(const TObject* iter);
+  Int_t RecV0(const TObject* iter);
+  Int_t RecBarrelTracks(const TObject* iter);
+
   /// pt cut for decay, minimum [GeV/c]
   float fPtMin;                                            //! transient
   /// Distance between decay tracks [cm] ??
@@ -111,9 +121,20 @@ class AliHLTD0Trigger : public AliHLTTrigger
   float fcosPoint;                                         //! transient 
 
   bool fplothisto;                                         //! transient 
+  bool fUseV0;                                             //! transient 
+
+  Double_t mD0PDG;                                         //! transient
 
   /// D0 inv. mass plot
-  TH1F *fD0mass;   //! transient  
+  TH1F *fD0mass;                                           //! transient  
+
+  vector<AliESDtrack*> fPos;                                //! transient
+  vector<AliESDtrack*> fNeg;                                //! transient
+
+  AliHLTD0toKpi *fd0calc;                                    //! transient
+  TObjArray *ftwoTrackArray;                                 //! transient
+
+  Int_t fTotalD0;                                            //! transient
 
   /// the default configuration entry for this component
   static const char* fgkOCDBEntry; //!transient
index 971c188..47705ec 100644 (file)
@@ -4,6 +4,11 @@
 #include "TMath.h"
 #include "AliESDtrack.h"
 #include "TVector3.h"
+#include "AliAODVertex.h"
+#include "AliESDVertex.h"
+#include "TObjArray.h"
+#include "AliVertexerTracks.h"
+
 
 ClassImp(AliHLTD0toKpi)
 
@@ -72,3 +77,38 @@ Double_t AliHLTD0toKpi::pointingAngle(AliESDtrack* n, AliESDtrack* p, Double_t *
 
   return TMath::Cos(pta); 
 }
+
+AliAODVertex* AliHLTD0toKpi::ReconstructSecondaryVertex(TObjArray *trkArray, Double_t b, const AliESDVertex *v)
+{
+  
+  AliESDVertex *vertexESD = 0;
+  AliAODVertex *vertexAOD = 0;
+  
+  AliVertexerTracks *vertexer = new AliVertexerTracks(b);
+  AliESDVertex* fV1 = new AliESDVertex(*v);
+  vertexer->SetVtxStart(fV1);
+  vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(trkArray);
+  delete vertexer; vertexer=NULL;
+  
+  if(!vertexESD) return vertexAOD;
+  
+  if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) { 
+    //AliDebug(2,"vertexing failed"); 
+    delete vertexESD; vertexESD=NULL;
+    return vertexAOD;
+  }
+
+  // convert to AliAODVertex
+  Double_t pos[3],cov[6],chi2perNDF;
+  vertexESD->GetXYZ(pos); // position
+  vertexESD->GetCovMatrix(cov); //covariance matrix
+  chi2perNDF = vertexESD->GetChi2toNDF();
+  //dispersion = vertexESD->GetDispersion();
+  delete vertexESD; vertexESD=NULL;
+
+  Int_t nprongs= trkArray->GetEntriesFast();
+  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
+
+  return vertexAOD;
+
+}
index 536f118..30035a3 100644 (file)
@@ -14,6 +14,9 @@
 #include "TObject.h"
 
 class AliESDtrack;
+class AliAODVertex;
+class AliESDVertex;
+class TObjArray;
 
 class AliHLTD0toKpi : public TObject
 {
@@ -22,6 +25,9 @@ public:
   Double_t InvMass(AliESDtrack* d1, AliESDtrack* d2);
   void cosThetaStar(AliESDtrack* n, AliESDtrack* p,Double_t &D0,Double_t &D0bar);
   Double_t pointingAngle(AliESDtrack* n, AliESDtrack* p, Double_t *pv, Double_t *sv);
+
+  AliAODVertex* ReconstructSecondaryVertex(TObjArray *trkArray, Double_t b, const AliESDVertex *v);
+
 private:
   
   ClassDef(AliHLTD0toKpi, 1)