]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/GammaConv/AliConversionTrackCuts.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliConversionTrackCuts.cxx
index a97c7ad427d8a1dce3d9f7de509c0c13b7770cfd..b911df725e1839815741cf2c97b44b5c33d52c28 100644 (file)
@@ -15,7 +15,7 @@
 
 /// @file   AliConversionTrackCuts.cxx
 /// @author Svein Lindal
-/// @brief  Base class for analysation of conversion particle - track correlations
+/// @brief  Base class for analysation of conversion particle - track correlations - track cuts
 
 
 #include "AliConversionTrackCuts.h"
@@ -25,6 +25,7 @@
 #include <iostream>
 #include "TH2F.h"
 #include "AliESDtrackCuts.h"
+#include "THn.h"
 
 using namespace std;
 ClassImp(AliConversionTrackCuts)
@@ -51,18 +52,21 @@ AliConversionTrackCuts::AliConversionTrackCuts() :
   fEsdTrackCutsExtra2(NULL),
   fEvent(NULL),
   fFilterBit(2048),
-  fDCAZmax(-1),
-  fDCAXYmax(-1),
+  fDCAZmax(3.2*3.2),
+  fDCAXYmax(2.4*2.4),
   fOwnedTracks(),
   fInitialized(kFALSE),
   fhPhi(NULL),
-  fhPt(NULL),
-  fhPhiPt(NULL),
+  //  fhPt(NULL),
+  //fhPhiPt(NULL),
   fhdcaxyPt(NULL),
   fhdcazPt(NULL),
   fhdca(NULL),
   fhnclpt(NULL),
   fhnclsfpt(NULL),
+  fhEtaPhi(NULL),
+  fhTrackEff(NULL),
+  fkCreateTrackEff(kFALSE),
   fHistograms(NULL) 
 {
   //Constructor
@@ -81,13 +85,16 @@ AliConversionTrackCuts::AliConversionTrackCuts(TString name, TString title = "ti
   fOwnedTracks(),
   fInitialized(kFALSE),
   fhPhi(NULL),  
-  fhPt(NULL),
-  fhPhiPt(NULL),
+  //fhPt(NULL),
+  //fhPhiPt(NULL),
   fhdcaxyPt(NULL),
   fhdcazPt(NULL),
   fhdca(NULL),
   fhnclpt(NULL),
   fhnclsfpt(NULL),
+  fhEtaPhi(NULL),
+  fhTrackEff(NULL),
+  fkCreateTrackEff(kFALSE),
   fHistograms(NULL)
 {
   //Constructor
@@ -125,25 +132,28 @@ void AliConversionTrackCuts::DefineESDCuts() {
   const Int_t filterbit = fFilterBit;
 
   if (filterbit == 128) {
-    fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
-    fEsdTrackCuts->SetMinNClustersTPC(70);
-    
+    if(!fEsdTrackCuts) {
+      fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+      fEsdTrackCuts->SetMinNClustersTPC(70);
+    }
   }  else if (filterbit == 256) {
-    // syst study
-    fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
-    fEsdTrackCuts->SetMinNClustersTPC(80);
-    fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
-    fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
-    fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
-    
+    if(!fEsdTrackCuts) {
+      // syst study
+      fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+      fEsdTrackCuts->SetMinNClustersTPC(80);
+      fEsdTrackCuts->SetMaxChi2PerClusterTPC(3);
+      fEsdTrackCuts->SetMaxDCAToVertexZ(2.7);
+      fEsdTrackCuts->SetMaxDCAToVertexXY(1.9);
+    }
   }  else if (filterbit == 512) {
-    // syst study
-    fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
-    fEsdTrackCuts->SetMinNClustersTPC(60);
-    fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
-    fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
-    fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
-    
+    if(!fEsdTrackCuts) {
+      // syst study
+      fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+      fEsdTrackCuts->SetMinNClustersTPC(60);
+      fEsdTrackCuts->SetMaxChi2PerClusterTPC(5);
+      fEsdTrackCuts->SetMaxDCAToVertexZ(3.7);
+      fEsdTrackCuts->SetMaxDCAToVertexXY(2.9);
+    }
   } else if (filterbit == 1024) {
     if(!fEsdTrackCuts) {
       fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
@@ -192,38 +202,74 @@ void AliConversionTrackCuts::DefineESDCuts() {
   }
 }
 
+
 //______________________________________________________________________________
 Bool_t AliConversionTrackCuts::AcceptTrack(AliESDtrack * track) {
   //Check esd track
+  FillHistograms(kPreCut, track);
+
+
+  if( fFilterBit == 256) {
+
+    ///Standalone tpc tracks constrained
+    const AliExternalTrackParam * param = track->GetConstrainedParam();
+    if(param) {
+      AliESDtrack* esdTrack = new AliESDtrack(*track);
+      esdTrack->CopyFromVTrack(param);
+      track = esdTrack;
+      fOwnedTracks.Add(track);
+
+      if( !fEsdTrackCuts->IsSelected(track)) return kFALSE;
+
+      FillHistograms(1, track);
+
+      Double_t dca[2];
+      GetDCA(track, dca);
+      
+      FillDCAHist(dca[1], dca[0], track);
+      if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
+      return kTRUE;
+    } else {
+      return kFALSE;
+    }
+
+    return kFALSE;
+  }
+
 
   if(!fInitialized) {
     DefineESDCuts();
-    if(fDCAXYmax > 0) {
-      if(fEsdTrackCuts) fEsdTrackCuts->SetMaxDCAToVertexXY(fDCAXYmax);
-    }
-    if(fDCAZmax > 0) {
-      if(fEsdTrackCuts) fEsdTrackCuts->SetMaxDCAToVertexZ(fDCAZmax);
-    }
+    // if(fDCAXYmax > 0) {
+    //   if(fEsdTrackCuts) fEsdTrackCuts->SetMaxDCAToVertexXY(fDCAXYmax);
+    // }
+    // if(fDCAZmax > 0) {
+    //   if(fEsdTrackCuts) fEsdTrackCuts->SetMaxDCAToVertexZ(fDCAZmax);
+    // }
   
     fInitialized = kTRUE;
   }
 
-  FillHistograms(kPreCut, track);
 
+  Double_t dca[2];
+  GetDCA(track, dca);
 
   
-  if( !fEsdTrackCuts->IsSelected(track)) return kFALSE;
-
-
   ///If only one track cuts then it has passed the cuts
   if( !(fEsdTrackCutsExtra1 && fEsdTrackCutsExtra2)) {
     FillHistograms(1, track);
+    FillDCAHist(dca[1], dca[0], track);
+    if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
     return kTRUE;
   }
+
   ///If passing extra
   if (fEsdTrackCutsExtra1 && fEsdTrackCutsExtra1->IsSelected(track)) {
-    FillHistograms(2, track);
     FillHistograms(1, track);
+    FillHistograms(2, track);
+
+    FillDCAHist(dca[1], dca[0], track);
+    if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
+    
     return kTRUE;
   } 
 
@@ -231,13 +277,17 @@ Bool_t AliConversionTrackCuts::AcceptTrack(AliESDtrack * track) {
   if (fEsdTrackCutsExtra2 && fEsdTrackCutsExtra2->IsSelected(track)) {
     const AliExternalTrackParam * param = track->GetConstrainedParam();
     if(param) {
-      AliESDtrack* esdTrack = new AliESDtrack(track);
+      AliESDtrack* esdTrack = new AliESDtrack(*track);
       esdTrack->CopyFromVTrack(param);
       track = esdTrack;
       fOwnedTracks.Add(track);
 
       FillHistograms(3, track);
       FillHistograms(1, track);
+
+      FillDCAHist(dca[1], dca[0], track);
+      if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
+
       return kTRUE;
     } else {
       return kFALSE;
@@ -259,70 +309,156 @@ Bool_t AliConversionTrackCuts::AcceptTrack(AliESDtrack * track) {
   // ///Get impact parameters
   // Double_t extCov[15];
   // track->GetExternalCovariance(extCov);
-  // Float_t b[2];
-  // Float_t bCov[3];
-  // track->GetImpactParameters(b,bCov);
-  // if (bCov[0]<=0 || bCov[2]<=0) {
-  //   AliDebug(1, "Estimated b resolution lower or equal zero!");
-  //   bCov[0]=0; bCov[2]=0;
-  // }
-  
-  // Float_t dcaToVertexXY = b[0];
-  // Float_t dcaToVertexZ = b[1];
-  // FillDCAHist(dcaToVertexZ, dcaToVertexXY, track);
   // return kTRUE;
 }
 
 Bool_t AliConversionTrackCuts::AcceptTrack(AliAODTrack * track) {
   //Check aod track
+  
   FillHistograms(kPreCut, track);
-  if(!track->TestFilterBit(fFilterBit)) {
-    return kFALSE;
-  }
+  
+  if (fFilterBit == 768) {
+    if(!track->IsHybridGlobalConstrainedGlobal()) return kFALSE;
+      
+    if (!(track->GetStatus() & AliVTrack::kITSrefit)) {
+      return kFALSE;
+    }
+      
+    //The cluster sharing cut can be done with:
+    Double_t frac = Double_t(track->GetTPCnclsS()) / Double_t(track->GetTPCncls());
+    if (frac > 0.4) return kFALSE;
+      
+    ///Do dca xy cut!
+    FillHistograms(1, track);
+      
+    ///DCA
+    Double_t dca[2] = { -999, -999};
+    //Bool_t dcaok = 
+    GetDCA(track, dca);
+    FillDCAHist(dca[1], dca[0], track);
+      
+      
+    if(track->IsGlobalConstrained()) {
+      FillHistograms(3, track);
+    } else {
+      FillHistograms(2, track);
+    }
+      
+    if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
+      
+    return kTRUE;
+    
+    ////////////////////////////////
+    //// Standalone
+    ////////////////////////////////
+  } else  if(fFilterBit == 256) {
+    if(!track->IsTPCConstrained()) return kFALSE;
 
-  ///Do dca xy cut!
 
-  return kTRUE;
 
+    ///DCA
+    Double_t dca[2] = { -999, -999};
+    GetDCA(track, dca);
 
-  // if (track->GetTPCNcls() < fTPCminNClusters) return kFALSE;
-  // FillHistograms(kCutNcls, track);
+    if( (dca[0]*dca[0]/fDCAXYmax + dca[1]*dca[1]/fDCAZmax) > 1 ) {
+      FillHistograms(3, track);
+      return kFALSE;
+    }
 
-  // if (track->Chi2perNDF() > fTPCmaxChi2) return kFALSE;
-  // FillHistograms(kCutNDF, track);
+    if(track->GetTPCncls() < 70) {
+      FillHistograms(4, track);
+      return kFALSE;
+    }
 
-  // AliAODVertex *vertex = track->GetProdVertex();
-  // if (vertex && fRejectKinkDaughters) {
-  //   if (vertex->GetType() == AliAODVertex::kKink) {
-  //     return kFALSE;
-  //   }
-  // }
-  // FillHistograms(kCutKinc, track);
-  
-  // if(TMath::Abs(track->ZAtDCA()) > fDCAZmax) {
-  //   return kFALSE;
-  // }
-  // FillHistograms(kCutDCAZ, track);
+    AliAODVertex * vtx = track->GetProdVertex();
+    if (vtx->GetType() == AliAODVertex::kKink ) {
+      FillHistograms(5, track);
+      return kFALSE;
+    }
 
-  // Float_t xatdca = track->XAtDCA();
-  // Float_t yatdca = track->YAtDCA();
-  // Float_t xy = xatdca*xatdca + yatdca*yatdca;
-  // if(xy > fDCAXYmax) {
-  //   return kFALSE;
-  // }
+    if(track->Chi2perNDF() > 36) {
+      FillHistograms(6, track);
+      return kFALSE;
+    }
+    if(track->Chi2perNDF() > 26) {
+      FillHistograms(7, track);
+      return kFALSE;
+    }
+    if(track->Chi2perNDF() > 16) {
+      FillHistograms(8, track);
+      return kFALSE;
+    }
+    if(track->Chi2perNDF() > 4) {
+      FillHistograms(9, track);
+      return kFALSE;
+    }
 
-  // FillHistograms(kCutDCAXY, track);
 
 
-  // fhnclpt->Fill(track->Pt(), track->GetTPCNcls());
-  // if(track->GetTPCNclsF() > 0) fhnclsfpt->Fill(track->Pt(), ((Double_t) track->GetTPCNcls())/track->GetTPCNclsF());
-  // FillDCAHist(track->ZAtDCA(), TMath::Sqrt(track->XAtDCA()*track->XAtDCA() + track->YAtDCA()*track->YAtDCA()), track);
+    FillDCAHist(dca[1], dca[0], track);
 
+    FillHistograms(2, track);
+    if(fhEtaPhi) fhEtaPhi->Fill(track->Eta(), track->Phi());
+    return kTRUE;
 
+  }
+  return kFALSE;
 }
 
 
+///______________________________________________________________________________
+Bool_t AliConversionTrackCuts::GetDCA(const AliESDtrack *track, Double_t dcaxyz[2]) {
+  ///Get track dca esd trck
+  Float_t dca[2];
+  Float_t bCov[3];
+  track->GetImpactParameters(dca,bCov);
+  if (bCov[0]<=0 || bCov[2]<=0) {
+    AliDebug(1, "Estimated b resolution lower or equal zero!");
+    bCov[0]=0; bCov[2]=0;
+    return kFALSE;
+  }
 
+  dcaxyz[0] = dca[0];
+  dcaxyz[1] = dca[1];
+  
+  return kTRUE;
+}
+
+///_____________________________________________________________________________
+Bool_t AliConversionTrackCuts::GetDCA(const AliAODTrack *track, Double_t dca[2]) {
+  ///Get track dca aod trck
+  if(track->TestBit(AliAODTrack::kIsDCA)){
+    dca[0]=track->DCA();
+    dca[1]=track->ZAtDCA();
+    return kTRUE;
+  }
+  
+  Bool_t ok=kFALSE;
+  if(fEvent) {
+    Double_t covdca[3];
+    //AliAODTrack copy(*track);
+    AliExternalTrackParam etp; etp.CopyFromVTrack(track);
+    
+    Float_t xstart = etp.GetX();
+    if(xstart>3.) {
+      dca[0]=-999.;
+      dca[1]=-999.;
+    //printf("This method can be used only for propagation inside the beam pipe \n");
+    return kFALSE;
+    }
+
+
+    AliAODVertex *vtx =(AliAODVertex*)(fEvent->GetPrimaryVertex());
+    Double_t fBzkG = fEvent->GetMagneticField(); // z componenent of field in kG
+    ok = etp.PropagateToDCA(vtx,fBzkG,kVeryBig,dca,covdca);
+    //ok = copy.PropagateToDCA(vtx,fBzkG,kVeryBig,dca,covdca);
+  }
+  if(!ok){
+    dca[0]=-999.;
+    dca[1]=-999.;
+  }
+  return ok;
+}
 
 
 
@@ -334,7 +470,7 @@ TList * AliConversionTrackCuts::CreateHistograms() {
   fHistograms->SetOwner(kTRUE);
   fHistograms->SetName("trackCuts");
 
-  fhPhi = new TH2F("phi", "phi", 5, -0.5, 4.5, 32, 0, TMath::TwoPi());
+  fhPhi = new TH2F(Form("phi_%s", GetName()), Form("phi_%s", GetTitle()), 5, -0.5, 4.5, 32, 0, TMath::TwoPi());
   // TAxis * xax = fhPhi->GetXaxis();
   // for(Int_t i = 0; i < kNCuts; i++){
   //   xax->SetBinLabel(xax->FindFixBin(i), fgkCutNames[i]);
@@ -342,6 +478,8 @@ TList * AliConversionTrackCuts::CreateHistograms() {
   fHistograms->Add(fhPhi);
   
 
+  fhEtaPhi = new TH2F(Form("etaphi_%s",GetName()), Form("etaphi_%s", GetTitle()), 36, -0.9, 0.9, 32, 0, TMath::TwoPi());
+  fHistograms->Add(fhEtaPhi);
 
   // fhPt = new TH2F("pt", "pt", kNCuts+2, kPreCut -0.5, kNCuts + 0.5, 
   //                             20, 0., 20.);
@@ -354,20 +492,41 @@ TList * AliConversionTrackCuts::CreateHistograms() {
   //  fhPhiPt = new TH2F("phipt", "phipt", 100, 0, 100, 64, 0, TMath::TwoPi());
   //fHistograms->Add(fhPhiPt);
 
-  // fhdcaxyPt = new TH2F("dcaxypt", "dcaxypt", 20, 0, 20, 50, 0, 5);
-  // fHistograms->Add(fhdcaxyPt);
+  fhdcaxyPt = new TH2F(Form("dcaxypt_%s", GetName()),  Form("dcaxypt_%s", GetTitle()), 20, 0, 20, 50, -2.5, 2.5);
+  fHistograms->Add(fhdcaxyPt);
+
+  fhdcazPt = new TH2F(Form("dcazpt_%s", GetName()),  Form("dcazpt_%s", GetTitle()), 20, 0, 20, 70, -3.5, 3.5);
+  fHistograms->Add(fhdcazPt);
 
-  // fhdcazPt = new TH2F("dcazpt", "dcazpt", 20, 0, 20, 50, 0, 5);
-  // fHistograms->Add(fhdcazPt);
+  fhdca = new TH2F(Form("dca_%s", GetName()),  Form("dca_%s", GetTitle()), 70, -3.5, 3.5, 50, -2.5, 2.5);
+  fhdca->SetXTitle("dca z");
+  fhdca->SetYTitle("dca xy");
 
-  // fhdca = new TH2F("dca", "dca", 60, -3, 3, 60, -3, 3);
-  // fHistograms->Add(fhdca);
+  
+  fHistograms->Add(fhdca);
 
   // fhnclpt = new TH2F("nclstpcvspt", "nclstpcvspt", 20, 0, 20, 50, 0, 100);
   // fHistograms->Add(fhnclpt);
 
   // fhnclsfpt = new TH2F("nclsfpt", "nclsfpt", 20, 0, 20, 60, 0, 1.2);
   // fHistograms->Add(fhnclsfpt);
+
+  
+
+  if (fkCreateTrackEff) {
+    const Double_t ptbins[23] = {0.5,  0.6, 0.7,  0.8, 0.9,  1.0, 1.2, 1.4,  1.6, 1.8, 
+                         2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.25, 
+                         4.5, 4.75, 5.0};
+    
+    const Int_t bins[4] = { 12,  22, 36, 32};
+    const Double_t min[4] = { -12, 0.5, -.9,  0};
+    const Double_t max[4] = {  12,   5,  .9,  2*TMath::Pi() };
+    
+    fhTrackEff = new THnF(Form("hTrackEff_%s", GetName()), "hTrackEff", 4, bins, min, max);
+    fhTrackEff->SetBinEdges(1, ptbins);
+    fHistograms->Add(fhTrackEff);
+  }
+
   
   return fHistograms;
 }
@@ -385,16 +544,21 @@ void AliConversionTrackCuts::FillDCAHist(Float_t dcaz, Float_t dcaxy, AliVTrack
   if(fhdcaxyPt) fhdcaxyPt->Fill(track->Pt(), dcaxy);
   if(fhdcazPt) fhdcazPt->Fill(track->Pt(), dcaz);
   if(fhdca) fhdca->Fill(dcaz, dcaxy);
+  
+  if(fhTrackEff) {
+    Double_t val[4];
+    val[0] = fEvent->GetPrimaryVertex()->GetZ();
+    val[1] =  track->Pt();
+    val[2] =  track->Eta();
+    val[3] =  track->Phi();
+    
+    fhTrackEff->Fill(val);
+  }
 }
 
 
 
 
-
-
-
-
-
 //_________________________________________________________________________________________________
 void AliConversionTrackCuts::Print(const Option_t *) const
 {