]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/vertexingHF/AliAODRecoCascadeHF.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAODRecoCascadeHF.cxx
index a4fd6fbd338abbd5a56c3709a9abe2826d7e28f2..ddfe39c1567ee5b8ba3e259594b87bfdadf939ca 100644 (file)
@@ -113,7 +113,7 @@ Double_t AliAODRecoCascadeHF::InvMassDstarKpipi() const
 //----------------------------------------------------------------------------
 Int_t AliAODRecoCascadeHF::MatchToMC(Int_t pdgabs,Int_t pdgabs2prong,
                                      Int_t *pdgDg,Int_t *pdgDg2prong,
-                                    TClonesArray *mcArray) const
+                                    TClonesArray *mcArray, Bool_t isV0) const
 {
   //
   // Check if this candidate is matched to a MC signal
@@ -122,32 +122,82 @@ Int_t AliAODRecoCascadeHF::MatchToMC(Int_t pdgabs,Int_t pdgabs2prong,
   // 
 
   Int_t ndg=GetNDaughters();
-  if(!ndg) {
+  if(ndg==0) {
     AliError("No daughters available");
     return -1;
   }
-  
-  AliAODRecoDecayHF2Prong *the2Prong = Get2Prong();
-  Int_t lab2Prong = the2Prong->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong);
+
+  if ( isV0 &&
+       ( (pdgDg[1]==2212 && pdgDg[0]==310) ||
+        (pdgDg[1]==211 && pdgDg[0]==3122) ) ) {
+    AliWarning("Please, pay attention: first element in AliAODRecoCascadeHF object must be the bachelor and second one V0. Skipping!");
+    return -1;
+  }
+
+  Int_t lab2Prong = -1;
+
+  if (!isV0) {
+    AliAODRecoDecayHF2Prong *the2Prong = Get2Prong();
+    lab2Prong = the2Prong->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong);
+  } else {
+    AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(Getv0());
+    lab2Prong = theV0->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong); // the V0
+  }
 
   if(lab2Prong<0) return -1;
 
   Int_t dgLabels[10]={0,0,0,0,0,0,0,0,0,0};
 
-  // loop on daughters and write labels
-  for(Int_t i=0; i<ndg; i++) {
-    AliVTrack *trk = (AliVTrack*)GetDaughter(i);
-    Int_t lab = trk->GetLabel();
-    if(lab==-1) { // this daughter is the 2prong
-      lab=lab2Prong;
-    } else if(lab<-1) {
-      printf("daughter with negative label\n");
-      continue;
+  if (!isV0) {
+    // loop on daughters and write labels
+    for(Int_t i=0; i<ndg; i++) {
+      AliVTrack *trk = dynamic_cast<AliVTrack*>(GetDaughter(i));
+      if(!trk) continue;
+      Int_t lab = trk->GetLabel();
+      if(lab==-1) { // this daughter is the 2prong
+       lab=lab2Prong;
+      } else if(lab<-1) continue;
+      dgLabels[i] = lab;
+    }
+  } else {
+    AliVTrack *trk = dynamic_cast<AliVTrack*>(GetBachelor()); // the bachelor
+    if (!trk) return -1;
+    dgLabels[0] = trk->GetLabel();//TMath::Abs(trk->GetLabel());
+    dgLabels[1] = lab2Prong;
+  }
+
+  Int_t finalLabel = AliAODRecoDecay::MatchToMC(pdgabs,mcArray,dgLabels,2,2,pdgDg);
+
+  if (finalLabel>=0){
+    // debug printouts for Lc->V0 bachelor case
+
+    if ( isV0 && (dgLabels[0]!=-1 && dgLabels[1]!=-1) ) {
+      AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(Getv0());
+      Bool_t onTheFly = theV0->GetOnFlyStatus();
+      if (pdgDg[0]==2212 && pdgDg[1]==310) {
+       AliAODMCParticle*k0s = dynamic_cast<AliAODMCParticle*>(mcArray->At(lab2Prong));
+       if(k0s){
+         Int_t labK0 = k0s->GetMother();       
+         AliAODMCParticle*k0bar = dynamic_cast<AliAODMCParticle*>(mcArray->At(labK0));
+         if(k0bar){
+           AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d) -> LabelK0S=%d (%d -> %d %d)",onTheFly,labK0,k0bar->GetPdgCode(),lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
+           AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
+                           finalLabel,pdgabs,
+                           dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
+         }
+       }
+      } else if (pdgDg[0]==211 && pdgDg[1]==3122) {
+       AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d -> %d %d)",onTheFly,lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
+       AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
+                       finalLabel,pdgabs,
+                     dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
+      }
+
     }
-    dgLabels[i] = lab;
   }
 
-  return AliAODRecoDecay::MatchToMC(pdgabs,mcArray,dgLabels,2,2,pdgDg);
+  return finalLabel;
+
 }
 //-----------------------------------------------------------------------------
 Bool_t AliAODRecoCascadeHF::SelectDstar(const Double_t *cutsDstar,
@@ -198,7 +248,7 @@ Bool_t AliAODRecoCascadeHF::SelectDstar(const Double_t *cutsDstar,
 }
 //-----------------------------------------------------------------------------
 Bool_t AliAODRecoCascadeHF::SelectLctoV0(const Double_t *cutsLctoV0, 
-                                        Bool_t okLck0sp, Bool_t okLcLpi) const 
+                                        Bool_t okLck0sp, Bool_t okLcLpi, Bool_t okLcLbarpi) const 
 {
   // cuts on Lambdac candidates to V0+bachelor
   // (to be passed to AliAODRecoDecayHF3Prong::SelectLctoV0())
@@ -209,14 +259,14 @@ Bool_t AliAODRecoCascadeHF::SelectLctoV0(const Double_t *cutsLctoV0,
   // 4 = pT min Bachelor track [GeV/c]
   // 5 = pT min V0-Positive track [GeV/c]
   // 6 = pT min V0-Negative track [GeV/c]
-  // 7 = dca cut on the V0 (cm)
-  // 8 = dca cut on the cascade (cm)
+  // 7 = dca cut on the cascade (cm)
+  // 8 = dca cut on the V0 (cm)
 
-//   if ( !Getv0() || !Getv0PositiveTrack() || !Getv0NegativeTrack() ) 
-//     { AliInfo(Form("Not adapted for ESDv0s, return true...")); return false; }
+  //   if ( !Getv0() || !Getv0PositiveTrack() || !Getv0NegativeTrack() ) 
+  //     { AliInfo(Form("Not adapted for ESDv0s, return true...")); return false; }
 
   Double_t mLck0sp,mLcLpi;
-  okLck0sp=1; okLcLpi=1;
+  okLck0sp=1; okLcLpi=1; okLcLbarpi=1;
   
   Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
   Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
@@ -236,26 +286,73 @@ Bool_t AliAODRecoCascadeHF::SelectLctoV0(const Double_t *cutsLctoV0,
   if(TMath::Abs(mLck0sp-mLcPDG)>cutsLctoV0[0]) okLck0sp = 0;
   //   with Lambda pi hypothesis
   if(TMath::Abs(mLcLpi-mLcPDG)>cutsLctoV0[1]) okLcLpi = 0;
-  
+  okLcLbarpi = okLcLpi;
   // cuts on the v0 mass
-  if(TMath::Abs(mk0s-mk0sPDG)>cutsLctoV0[2]) okLck0sp = 0;
-  if( TMath::Abs(mlambda-mLPDG)>cutsLctoV0[3] && 
-      TMath::Abs(malambda-mLPDG)>cutsLctoV0[3] ) okLcLpi = 0;
+  if( TMath::Abs(mk0s-mk0sPDG)>cutsLctoV0[2]) okLck0sp = 0;
+  //if( TMath::Abs(mlambda-mLPDG)>cutsLctoV0[3] && 
+  //TMath::Abs(malambda-mLPDG)>cutsLctoV0[3] ) okLcLpi = 0;
+  if( !(GetBachelor()->Charge()==+1 && TMath::Abs(mlambda-mLPDG)<=cutsLctoV0[3]) ) okLcLpi = 0;
+  if( !(GetBachelor()->Charge()==-1 && TMath::Abs(malambda-mLPDG)<=cutsLctoV0[3]) ) okLcLbarpi = 0;
   
-  if(!okLck0sp && !okLcLpi) return 0;
+  if(!okLck0sp && !okLcLpi && !okLcLbarpi) return 0;
   
   // cuts on the minimum pt of the tracks 
   if(TMath::Abs(GetBachelor()->Pt()) < cutsLctoV0[4]) return 0;
   if(TMath::Abs(Getv0PositiveTrack()->Pt()) < cutsLctoV0[5]) return 0;
   if(TMath::Abs(Getv0NegativeTrack()->Pt()) < cutsLctoV0[6]) return 0;
   
-  // cut on the v0 dca
-  if(TMath::Abs(Getv0()->DcaV0Daughters()) > cutsLctoV0[7]) return 0;
-  
   // cut on the cascade dca
-  if( TMath::Abs(GetDCA(0))>cutsLctoV0[8] ||
-      TMath::Abs(Getv0()->DcaPosToPrimVertex())>cutsLctoV0[8] ||
-      TMath::Abs(Getv0()->DcaNegToPrimVertex())>cutsLctoV0[8] ) return 0;
+  if( TMath::Abs(GetDCA(0))>cutsLctoV0[7] //||
+      //TMath::Abs(Getv0()->DcaPosToPrimVertex())>cutsLctoV0[7] ||
+      //TMath::Abs(Getv0()->DcaNegToPrimVertex())>cutsLctoV0[7]
+      ) return 0;
+  
+  // cut on the v0 dca
+  if(TMath::Abs(Getv0()->DcaV0Daughters()) > cutsLctoV0[8]) return 0;
+
+  // cut on V0 cosine of pointing angle wrt PV
+  if (CosV0PointingAngle() < cutsLctoV0[9]) { // cosine of V0 pointing angle wrt primary vertex
+    AliDebug(4,Form(" V0 cosine of pointing angle doesn't pass the cut"));
+    return 0;
+  }
+
+  // cut on bachelor transverse impact parameter wrt PV
+  if (TMath::Abs(Getd0Prong(0)) > cutsLctoV0[10]) { // bachelor transverse impact parameter wrt PV
+    AliDebug(4,Form(" bachelor transverse impact parameter doesn't pass the cut"));
+    return 0;
+  }
+
+  // cut on V0 transverse impact parameter wrt PV
+  if (TMath::Abs(Getd0Prong(1)) > cutsLctoV0[11]) { // V0 transverse impact parameter wrt PV
+    AliDebug(4,Form(" V0 transverse impact parameter doesn't pass the cut"));
+    return 0;
+  }
+
+  // cut on K0S invariant mass veto
+  if (TMath::Abs(Getv0()->MassK0Short()-mk0sPDG) < cutsLctoV0[12]) { // K0S invariant mass veto
+    AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
+    return 0;
+  }
+
+  // cut on Lambda/LambdaBar invariant mass veto
+  if (TMath::Abs(Getv0()->MassLambda()-mLPDG) < cutsLctoV0[13] ||
+      TMath::Abs(Getv0()->MassAntiLambda()-mLPDG) < cutsLctoV0[13] ) { // Lambda/LambdaBar invariant mass veto
+    AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
+    return 0;
+  }
+
+  // cut on gamma invariant mass veto                                                                                                                      
+  if (Getv0()->InvMass2Prongs(0,1,11,11) < cutsLctoV0[14]) { // K0S invariant mass veto
+    AliDebug(4,Form(" veto on gamma invariant mass doesn't pass the cut"));
+    return 0;
+  }
+
+  // cut on V0 pT min                                                                                                                                      
+  if (Getv0()->Pt() < cutsLctoV0[15]) { // V0 pT min                                                                                         
+    AliDebug(4,Form(" V0 track Pt=%2.2e > %2.2e",Getv0()->Pt(),cutsLctoV0[15]));
+    return 0;
+  }
   
   return true; 
 
@@ -299,3 +396,108 @@ Bool_t AliAODRecoCascadeHF::TrigonometricalCut() const {
   if((phi01>phi00) || (phi02>phi00)) return kFALSE;
   return kTRUE;
 }
+
+//-----------------------------------------------------------------------------
+Double_t AliAODRecoCascadeHF::DecayLengthV0() const
+{
+  //
+  // Returns V0 decay length wrt primary vertex
+  //
+
+  AliAODv0 *v0 = (AliAODv0*)Getv0();
+
+  if (!v0) 
+    return -1.;
+  AliAODVertex *vtxPrimary = GetPrimaryVtx();
+  Double_t posVtx[3] = {0.,0.,0.};
+  vtxPrimary->GetXYZ(posVtx);
+  return v0->DecayLengthV0(posVtx);
+
+}
+//-----------------------------------------------------------------------------
+Double_t AliAODRecoCascadeHF::DecayLengthXYV0() const
+{
+  //
+  // Returns transverse V0 decay length wrt primary vertex
+  //
+  AliAODv0 *v0 = (AliAODv0*)Getv0();
+
+  if (!v0) 
+    return -1.;
+  AliAODVertex *vtxPrimary = GetPrimaryVtx();
+  Double_t posVtx[3] = {0.,0.,0.};
+  vtxPrimary->GetXYZ(posVtx);
+  return v0->DecayLengthXY(posVtx);
+
+}
+//-----------------------------------------------------------------------------
+Double_t AliAODRecoCascadeHF::CosV0PointingAngle() const 
+{
+  //
+  // Returns cosine of V0 pointing angle wrt primary vertex
+  //
+
+  AliAODv0 *v0 = (AliAODv0*)Getv0();
+
+  if (!v0) 
+    return -999.;
+
+  AliAODVertex *vtxPrimary = GetPrimaryVtx();
+  Double_t posVtx[3] = {0.,0.,0.};
+  vtxPrimary->GetXYZ(posVtx);
+  return v0->CosPointingAngle(posVtx);
+
+}
+//-----------------------------------------------------------------------------
+Double_t AliAODRecoCascadeHF::CosV0PointingAngleXY() const 
+{
+  //
+  // Returns XY cosine of V0 pointing angle wrt primary vertex
+  //
+
+  AliAODv0 *v0 = (AliAODv0*)Getv0();
+
+  if (!v0) 
+    return -999.;
+
+  AliAODVertex *vtxPrimary = GetPrimaryVtx();
+  Double_t posVtx[3] = {0.,0.,0.};
+  vtxPrimary->GetXYZ(posVtx);
+  return v0->CosPointingAngleXY(posVtx);
+
+}
+//-----------------------------------------------------------------------------
+Double_t AliAODRecoCascadeHF::NormalizedV0DecayLength() const
+{
+  //
+  // Returns V0 normalized decay length wrt primary vertex
+  //
+
+  AliAODv0 *v0 = (AliAODv0*)Getv0();
+
+  if (!v0) 
+    return -1.;
+  //AliAODVertex *vtxPrimary = GetPrimaryVtx();
+  //Double_t posVtx[3] = {0.,0.,0.};
+  //vtxPrimary->GetXYZ(posVtx);
+  //return v0->NormalizedDecayLength(posVtx);
+  return v0->NormalizedDecayLength(GetPrimaryVtx());
+
+}
+//-----------------------------------------------------------------------------
+Double_t AliAODRecoCascadeHF::NormalizedV0DecayLengthXY() const
+{
+  //
+  // Returns transverse V0 normalized decay length wrt primary vertex
+  //
+  AliAODv0 *v0 = (AliAODv0*)Getv0();
+
+  if (!v0) 
+    return -1.;
+  //AliAODVertex *vtxPrimary = GetPrimaryVtx();
+  //Double_t posVtx[3] = {0.,0.,0.};
+  //vtxPrimary->GetXYZ(posVtx);
+  //return v0->NormalizedDecayLengthXY(posVtx);
+  return v0->NormalizedDecayLengthXY(GetPrimaryVtx());
+
+}