Changes needed to properly add and subtract the SDD low threshold with data compresse...

[u/mrichter/AliRoot.git] / PWG1 / AliComparisonDCA.cxx

//------------------------------------------------------------------------------\r
// Implementation of AliComparisonDCA class. It keeps information from \r
// comparison of reconstructed and MC particle tracks. In addtion, \r
// it keeps selection cuts used during comparison. The comparison \r
// information is stored in the ROOT histograms. Analysis of these \r
// histograms can be done by using Analyse() class function. The result of \r
// the analysis (histograms) are stored in the output picture_dca.root file.\r
//  \r
// Author: J.Otwinowski 04/02/2008 \r
//------------------------------------------------------------------------------\r
\r
/*\r
  // after running analysis, read the file, and get component\r
  gSystem->Load("libPWG1.so");\r
  TFile f("Output.root");\r
  AliComparisonDCA * comp = (AliComparisonDCA*)f.Get("AliComparisonDCA");\r
\r
  // analyse comparison data (output stored in pictures_dca.root)\r
  comp->Analyse();\r
\r
  // TPC track length parameterisation\r
  TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // length function\r
  TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);\r
  fl2.SetParameter(1,1);\r
  fl2.SetParameter(0,1);\r
*/\r
\r
#include <iostream>\r
\r
#include "TFile.h"\r
#include "TCint.h"\r
#include "TH3F.h"\r
#include "TH2F.h"\r
#include "TF1.h"\r
#include "TProfile.h"\r
#include "TProfile2D.h"\r
#include "TGraph2D.h"\r
#include "TCanvas.h"\r
#include "TGraph.h"\r
// \r
#include "AliTracker.h"   \r
#include "AliESDEvent.h"   \r
#include "AliESD.h"\r
#include "AliESDfriend.h"\r
#include "AliESDfriendTrack.h"\r
#include "AliRecInfoCuts.h" \r
#include "AliMCInfoCuts.h" \r
#include "AliLog.h" \r
#include "AliESDVertex.h" \r
//\r
#include "AliMathBase.h"\r
#include "AliTreeDraw.h" \r
\r
#include "AliMCInfo.h" \r
#include "AliESDRecInfo.h" \r
#include "AliComparisonDCA.h" \r
\r
using namespace std;\r
\r
ClassImp(AliComparisonDCA)\r
\r
//_____________________________________________________________________________\r
AliComparisonDCA::AliComparisonDCA():\r
  TNamed("AliComparisonDCA","AliComparisonDCA"),\r
\r
  // DCA histograms\r
  fD0TanSPtB1(0),\r
  fD1TanSPtB1(0),\r
  fD0TanSPtL1(0),\r
  fD1TanSPtL1(0),\r
  fD0TanSPtInTPC(0),\r
  fD1TanSPtInTPC(0),\r
  fVertex(0),\r
\r
  // Cuts \r
  fCutsRC(0), \r
  fCutsMC(0)  \r
{\r
  InitHisto();\r
  InitCuts();\r
\r
  // vertex (0,0,0)\r
  fVertex = new AliESDVertex();\r
  fVertex->SetXv(0.0);\r
  fVertex->SetYv(0.0);\r
  fVertex->SetZv(0.0);\r
}\r
\r
//_____________________________________________________________________________\r
AliComparisonDCA::~AliComparisonDCA()\r
{\r
  //\r
  if(fD0TanSPtB1) delete fD0TanSPtB1; fD0TanSPtB1=0;\r
  if(fD1TanSPtB1) delete fD1TanSPtB1; fD1TanSPtB1=0;\r
  if(fD0TanSPtL1) delete fD0TanSPtL1; fD0TanSPtL1=0;\r
  if(fD1TanSPtL1) delete fD1TanSPtL1; fD1TanSPtL1=0;\r
  if(fD0TanSPtInTPC) delete fD0TanSPtInTPC; fD0TanSPtInTPC=0;\r
  if(fD1TanSPtInTPC) delete fD1TanSPtInTPC; fD1TanSPtInTPC=0;\r
  if(fVertex) delete fVertex; fVertex=0;\r
}\r
\r
//_____________________________________________________________________________\r
void AliComparisonDCA::InitHisto()\r
{\r
  // DCA histograms\r
  fD0TanSPtB1 = new TH3F("DCAyTanSPtB1","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);\r
  fD0TanSPtB1->SetXTitle("tan(#theta)");\r
  fD0TanSPtB1->SetYTitle("#sqrt{p_{t}(GeV/c)}");\r
  fD0TanSPtB1->SetZTitle("DCA_{xy}");\r
\r
  fD1TanSPtB1 = new TH3F("DCAzTanSPtB1","DCAzTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);\r
  fD1TanSPtB1->SetXTitle("tan(#theta)");\r
  fD1TanSPtB1->SetYTitle("#sqrt(p_{t}(GeV/c))");\r
  fD1TanSPtB1->SetZTitle("DCA_{z}");\r
\r
  fD0TanSPtL1 = new TH3F("DCAyTanSPtL1","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);\r
  fD0TanSPtL1->SetXTitle("tan(#theta)");\r
  fD0TanSPtL1->SetYTitle("#sqrt{p_{t}(GeV/c)}");\r
  fD0TanSPtL1->SetZTitle("DCA_{xy}");\r
\r
  fD1TanSPtL1 = new TH3F("DCAzTanSPtL1","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1);\r
  fD1TanSPtL1->SetXTitle("tan(#theta)");\r
  fD1TanSPtL1->SetYTitle("#sqrt{p_{t}(GeV/c)}");\r
  fD1TanSPtL1->SetZTitle("DCA_{z}");\r
\r
  fD0TanSPtInTPC = new TH3F("DCAyTanSPtInTPC","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);\r
  fD0TanSPtInTPC->SetXTitle("tan(#theta)");\r
  fD0TanSPtInTPC->SetYTitle("#sqrt{p_{t}(GeV/c)}");\r
  fD0TanSPtInTPC->SetZTitle("DCA_{xy}");\r
\r
  fD1TanSPtInTPC = new TH3F("DCAzTanSPtInTPC","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1);\r
  fD1TanSPtInTPC->SetXTitle("tan(#theta)");\r
  fD1TanSPtInTPC->SetYTitle("#sqrt{p_{t}(GeV/c)}");\r
  fD1TanSPtInTPC->SetZTitle("DCA_{z}");\r
}\r
\r
//_____________________________________________________________________________\r
void AliComparisonDCA::InitCuts()\r
{\r
  if(!fCutsMC) \r
    AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");\r
  if(!fCutsRC) \r
    AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");\r
}\r
\r
//_____________________________________________________________________________\r
void AliComparisonDCA::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)\r
{\r
  // Fill DCA comparison information\r
  AliExternalTrackParam *track = 0;\r
  Double_t kRadius    = 3.0;      // beam pipe radius\r
  Double_t kMaxStep   = 5.0;      // max step\r
  Double_t field      = AliTracker::GetBz(); // nominal Bz field [kG]\r
  Double_t kMaxD      = 123456.0; // max distance\r
\r
  Int_t clusterITS[200];\r
  Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z\r
\r
  Float_t mcpt = infoMC->GetParticle().Pt();\r
  Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);\r
  Float_t spt = TMath::Sqrt(mcpt);\r
\r
  // distance to Prim. vertex \r
  const Double_t* dv = infoMC->GetVDist(); \r
\r
  Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();\r
\r
  // Check selection cuts\r
  if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; \r
  if (!isPrim) return;\r
  if (infoRC->GetStatus(1)!=3) return;\r
  if (!infoRC->GetESDtrack()) return;  \r
  if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;\r
  if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;\r
\r
  // calculate and set prim. vertex\r
  fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );\r
  fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );\r
  fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );\r
\r
  // calculate track parameters at vertex\r
  if (infoRC->GetESDtrack()->GetTPCInnerParam())\r
  {\r
    if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )\r
    {\r
      Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);\r
      Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);\r
\r
      if(bStatus && bDCAStatus) {\r
        fD0TanSPtInTPC->Fill(tantheta,spt,dca[0]);\r
        fD1TanSPtInTPC->Fill(tantheta,spt,dca[1]);\r
	  }\r
\r
	  delete track;\r
    }\r
  }\r
\r
 if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)==0){\r
    fD0TanSPtB1->Fill(tantheta,spt,dca[0]);\r
    fD1TanSPtB1->Fill(tantheta,spt,dca[1]);\r
  }\r
    fD0TanSPtL1->Fill(tantheta,spt,dca[0]);\r
    fD1TanSPtL1->Fill(tantheta,spt,dca[1]);\r
}\r
\r
//_____________________________________________________________________________\r
Long64_t AliComparisonDCA::Merge(TCollection* list) \r
{\r
  // Merge list of objects (needed by PROOF)\r
\r
  if (!list)\r
  return 0;\r
\r
  if (list->IsEmpty())\r
  return 1;\r
\r
  TIterator* iter = list->MakeIterator();\r
  TObject* obj = 0;\r
\r
  // collection of generated histograms\r
  Int_t count=0;\r
  while((obj = iter->Next()) != 0) \r
  {\r
    AliComparisonDCA* entry = dynamic_cast<AliComparisonDCA*>(obj);\r
    if (entry == 0) { \r
      Error("Add","Attempt to add object of class: %s to a %s",\r
	  obj->ClassName(),this->ClassName());\r
      return -1;\r
    }\r
\r
    fD0TanSPtB1->Add(entry->fD0TanSPtB1);\r
    fD1TanSPtB1->Add(entry->fD1TanSPtB1);\r
    fD0TanSPtL1->Add(entry->fD0TanSPtL1);\r
    fD1TanSPtL1->Add(entry->fD1TanSPtL1);\r
    fD0TanSPtInTPC->Add(entry->fD0TanSPtInTPC);\r
    fD1TanSPtInTPC->Add(entry->fD1TanSPtInTPC);\r
\r
    count++;\r
  }\r
\r
return count;\r
}\r
\r
//_____________________________________________________________________________\r
void AliComparisonDCA::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){\r
  // Process comparison information\r
  Process(infoMC,infoRC);\r
}\r
\r
//_____________________________________________________________________________\r
void AliComparisonDCA::Analyse()\r
{\r
  // Analyse output histograms\r
  \r
  AliComparisonDCA * comp=this;\r
  TGraph2D * gr=0;\r
  TGraph * gr0=0;\r
\r
  TFile *fp = new TFile("pictures_dca.root","recreate");\r
  fp->cd();\r
\r
  TCanvas * c = new TCanvas("DCA","DCA resloution");\r
  c->cd();\r
\r
  // DCA resolution\r
  gr0 = AliMathBase::MakeStat1D(comp->fD0TanSPtB1,2,5);\r
  gr0->GetXaxis()->SetTitle("Tan(#theta)");\r
  gr0->GetYaxis()->SetTitle("#sigmaDCA (cm)");\r
  gr0->Write("DCAResolTan");\r
  //\r
  gr = AliMathBase::MakeStat2D(comp->fD0TanSPtB1,4,2,5); \r
  gr->GetXaxis()->SetTitle("Tan(#theta)");\r
  gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV/c)}");\r
  gr->GetYaxis()->SetTitle("#sigmaDCA (cm)");\r
  gr->GetHistogram()->Write("DCAResolSPTTan");\r
\r
  gPad->Clear();\r
  gr0->Draw("al*");\r
\r
  fp->Close();\r
}\r