PWG1/AliComparisonDCA.cxx

   1 //------------------------------------------------------------------------------\r
   2 // Implementation of AliComparisonDCA class. It keeps information from \r
   3 // comparison of reconstructed and MC particle tracks. In addtion, \r
   4 // it keeps selection cuts used during comparison. The comparison \r
   5 // information is stored in the ROOT histograms. Analysis of these \r
   6 // histograms can be done by using Analyse() class function. The result of \r
   7 // the analysis (histograms) are stored in the output picture_dca.root file.\r
   8 //  \r
   9 // Author: J.Otwinowski 04/02/2008 \r
  10 //------------------------------------------------------------------------------\r
  11 \r
  12 /*\r
  13   // after running analysis, read the file, and get component\r
  14   gSystem->Load("libPWG1.so");\r
  15   TFile f("Output.root");\r
  16   AliComparisonDCA * comp = (AliComparisonDCA*)f.Get("AliComparisonDCA");\r
  17 \r
  18   // analyse comparison data (output stored in pictures_dca.root)\r
  19   comp->Analyse();\r
  20 \r
  21   // TPC track length parameterisation\r
  22   TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // length function\r
  23   TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);\r
  24   fl2.SetParameter(1,1);\r
  25   fl2.SetParameter(0,1);\r
  26 */\r
  27 \r
  28 #include <iostream>\r
  29 \r
  30 #include "TFile.h"\r
  31 #include "TCint.h"\r
  32 #include "TH3F.h"\r
  33 #include "TH2F.h"\r
  34 #include "TF1.h"\r
  35 #include "TProfile.h"\r
  36 #include "TProfile2D.h"\r
  37 #include "TGraph2D.h"\r
  38 #include "TCanvas.h"\r
  39 #include "TGraph.h"\r
  40 // \r
  41 #include "AliTracker.h"   \r
  42 #include "AliESDEvent.h"   \r
  43 #include "AliESD.h"\r
  44 #include "AliESDfriend.h"\r
  45 #include "AliESDfriendTrack.h"\r
  46 #include "AliRecInfoCuts.h" \r
  47 #include "AliMCInfoCuts.h" \r
  48 #include "AliLog.h" \r
  49 #include "AliESDVertex.h" \r
  50 //\r
  51 #include "AliMathBase.h"\r
  52 #include "AliTreeDraw.h" \r
  53 \r
  54 #include "AliMCInfo.h" \r
  55 #include "AliESDRecInfo.h" \r
  56 #include "AliComparisonDCA.h" \r
  57 \r
  58 using namespace std;\r
  59 \r
  60 ClassImp(AliComparisonDCA)\r
  61 \r
  62 //_____________________________________________________________________________\r
  63 AliComparisonDCA::AliComparisonDCA():\r
  64   TNamed("AliComparisonDCA","AliComparisonDCA"),\r
  65 \r
  66   // DCA histograms\r
  67   fD0TanSPtB1(0),\r
  68   fD1TanSPtB1(0),\r
  69   fD0TanSPtL1(0),\r
  70   fD1TanSPtL1(0),\r
  71   fD0TanSPtInTPC(0),\r
  72   fD1TanSPtInTPC(0),\r
  73   fVertex(0),\r
  74 \r
  75   // Cuts \r
  76   fCutsRC(0), \r
  77   fCutsMC(0)  \r
  78 {\r
  79   InitHisto();\r
  80   InitCuts();\r
  81 \r
  82   // vertex (0,0,0)\r
  83   fVertex = new AliESDVertex();\r
  84   fVertex->SetXv(0.0);\r
  85   fVertex->SetYv(0.0);\r
  86   fVertex->SetZv(0.0);\r
  87 }\r
  88 \r
  89 //_____________________________________________________________________________\r
  90 AliComparisonDCA::~AliComparisonDCA()\r
  91 {\r
  92   //\r
  93   if(fD0TanSPtB1) delete fD0TanSPtB1; fD0TanSPtB1=0;\r
  94   if(fD1TanSPtB1) delete fD1TanSPtB1; fD1TanSPtB1=0;\r
  95   if(fD0TanSPtL1) delete fD0TanSPtL1; fD0TanSPtL1=0;\r
  96   if(fD1TanSPtL1) delete fD1TanSPtL1; fD1TanSPtL1=0;\r
  97   if(fD0TanSPtInTPC) delete fD0TanSPtInTPC; fD0TanSPtInTPC=0;\r
  98   if(fD1TanSPtInTPC) delete fD1TanSPtInTPC; fD1TanSPtInTPC=0;\r
  99   if(fVertex) delete fVertex; fVertex=0;\r
 100 }\r
 101 \r
 102 //_____________________________________________________________________________\r
 103 void AliComparisonDCA::InitHisto()\r
 104 {\r
 105   // DCA histograms\r
 106   fD0TanSPtB1 = new TH3F("DCAyTanSPtB1","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);\r
 107   fD0TanSPtB1->SetXTitle("tan(#theta)");\r
 108   fD0TanSPtB1->SetYTitle("#sqrt{p_{t}(GeV/c)}");\r
 109   fD0TanSPtB1->SetZTitle("DCA_{xy}");\r
 110 \r
 111   fD1TanSPtB1 = new TH3F("DCAzTanSPtB1","DCAzTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);\r
 112   fD1TanSPtB1->SetXTitle("tan(#theta)");\r
 113   fD1TanSPtB1->SetYTitle("#sqrt(p_{t}(GeV/c))");\r
 114   fD1TanSPtB1->SetZTitle("DCA_{z}");\r
 115 \r
 116   fD0TanSPtL1 = new TH3F("DCAyTanSPtL1","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);\r
 117   fD0TanSPtL1->SetXTitle("tan(#theta)");\r
 118   fD0TanSPtL1->SetYTitle("#sqrt{p_{t}(GeV/c)}");\r
 119   fD0TanSPtL1->SetZTitle("DCA_{xy}");\r
 120 \r
 121   fD1TanSPtL1 = new TH3F("DCAzTanSPtL1","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1);\r
 122   fD1TanSPtL1->SetXTitle("tan(#theta)");\r
 123   fD1TanSPtL1->SetYTitle("#sqrt{p_{t}(GeV/c)}");\r
 124   fD1TanSPtL1->SetZTitle("DCA_{z}");\r
 125 \r
 126   fD0TanSPtInTPC = new TH3F("DCAyTanSPtInTPC","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);\r
 127   fD0TanSPtInTPC->SetXTitle("tan(#theta)");\r
 128   fD0TanSPtInTPC->SetYTitle("#sqrt{p_{t}(GeV/c)}");\r
 129   fD0TanSPtInTPC->SetZTitle("DCA_{xy}");\r
 130 \r
 131   fD1TanSPtInTPC = new TH3F("DCAzTanSPtInTPC","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1);\r
 132   fD1TanSPtInTPC->SetXTitle("tan(#theta)");\r
 133   fD1TanSPtInTPC->SetYTitle("#sqrt{p_{t}(GeV/c)}");\r
 134   fD1TanSPtInTPC->SetZTitle("DCA_{z}");\r
 135 }\r
 136 \r
 137 //_____________________________________________________________________________\r
 138 void AliComparisonDCA::InitCuts()\r
 139 {\r
 140   if(!fCutsMC) \r
 141     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");\r
 142   if(!fCutsRC) \r
 143     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");\r
 144 }\r
 145 \r
 146 //_____________________________________________________________________________\r
 147 void AliComparisonDCA::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)\r
 148 {\r
 149   // Fill DCA comparison information\r
 150   AliExternalTrackParam *track = 0;\r
 151   Double_t kRadius    = 3.0;      // beam pipe radius\r
 152   Double_t kMaxStep   = 5.0;      // max step\r
 153   Double_t field      = AliTracker::GetBz(); // nominal Bz field [kG]\r
 154   Double_t kMaxD      = 123456.0; // max distance\r
 155 \r
 156   Int_t clusterITS[200];\r
 157   Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z\r
 158 \r
 159   Float_t mcpt = infoMC->GetParticle().Pt();\r
 160   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);\r
 161   Float_t spt = TMath::Sqrt(mcpt);\r
 162 \r
 163   // distance to Prim. vertex \r
 164   const Double_t* dv = infoMC->GetVDist(); \r
 165 \r
 166   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();\r
 167 \r
 168   // Check selection cuts\r
 169   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; \r
 170   if (!isPrim) return;\r
 171   if (infoRC->GetStatus(1)!=3) return;\r
 172   if (!infoRC->GetESDtrack()) return;  \r
 173   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;\r
 174   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;\r
 175 \r
 176   // calculate and set prim. vertex\r
 177   fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );\r
 178   fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );\r
 179   fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );\r
 180 \r
 181   // calculate track parameters at vertex\r
 182   if (infoRC->GetESDtrack()->GetTPCInnerParam())\r
 183   {\r
 184     if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )\r
 185     {\r
 186       Bool_t bStatus = AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);\r
 187       Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);\r
 188 \r
 189       if(bStatus && bDCAStatus) {\r
 190         fD0TanSPtInTPC->Fill(tantheta,spt,dca[0]);\r
 191         fD1TanSPtInTPC->Fill(tantheta,spt,dca[1]);\r
 192           }\r
 193 \r
 194           delete track;\r
 195     }\r
 196   }\r
 197 \r
 198  if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)==0){\r
 199     fD0TanSPtB1->Fill(tantheta,spt,dca[0]);\r
 200     fD1TanSPtB1->Fill(tantheta,spt,dca[1]);\r
 201   }\r
 202     fD0TanSPtL1->Fill(tantheta,spt,dca[0]);\r
 203     fD1TanSPtL1->Fill(tantheta,spt,dca[1]);\r
 204 }\r
 205 \r
 206 //_____________________________________________________________________________\r
 207 Long64_t AliComparisonDCA::Merge(TCollection* list) \r
 208 {\r
 209   // Merge list of objects (needed by PROOF)\r
 210 \r
 211   if (!list)\r
 212   return 0;\r
 213 \r
 214   if (list->IsEmpty())\r
 215   return 1;\r
 216 \r
 217   TIterator* iter = list->MakeIterator();\r
 218   TObject* obj = 0;\r
 219 \r
 220   // collection of generated histograms\r
 221   Int_t count=0;\r
 222   while((obj = iter->Next()) != 0) \r
 223   {\r
 224     AliComparisonDCA* entry = dynamic_cast<AliComparisonDCA*>(obj);\r
 225     if (entry == 0) { \r
 226       Error("Add","Attempt to add object of class: %s to a %s",\r
 227           obj->ClassName(),this->ClassName());\r
 228       return -1;\r
 229     }\r
 230 \r
 231     fD0TanSPtB1->Add(entry->fD0TanSPtB1);\r
 232     fD1TanSPtB1->Add(entry->fD1TanSPtB1);\r
 233     fD0TanSPtL1->Add(entry->fD0TanSPtL1);\r
 234     fD1TanSPtL1->Add(entry->fD1TanSPtL1);\r
 235     fD0TanSPtInTPC->Add(entry->fD0TanSPtInTPC);\r
 236     fD1TanSPtInTPC->Add(entry->fD1TanSPtInTPC);\r
 237 \r
 238     count++;\r
 239   }\r
 240 \r
 241 return count;\r
 242 }\r
 243 \r
 244 //_____________________________________________________________________________\r
 245 void AliComparisonDCA::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){\r
 246   // Process comparison information\r
 247   Process(infoMC,infoRC);\r
 248 }\r
 249 \r
 250 //_____________________________________________________________________________\r
 251 void AliComparisonDCA::Analyse()\r
 252 {\r
 253   // Analyse output histograms\r
 254   \r
 255   AliComparisonDCA * comp=this;\r
 256   TGraph2D * gr=0;\r
 257   TGraph * gr0=0;\r
 258 \r
 259   TFile *fp = new TFile("pictures_dca.root","recreate");\r
 260   fp->cd();\r
 261 \r
 262   TCanvas * c = new TCanvas("DCA","DCA resloution");\r
 263   c->cd();\r
 264 \r
 265   // DCA resolution\r
 266   gr0 = AliMathBase::MakeStat1D(comp->fD0TanSPtB1,2,5);\r
 267   gr0->GetXaxis()->SetTitle("Tan(#theta)");\r
 268   gr0->GetYaxis()->SetTitle("#sigmaDCA (cm)");\r
 269   gr0->Write("DCAResolTan");\r
 270   //\r
 271   gr = AliMathBase::MakeStat2D(comp->fD0TanSPtB1,4,2,5); \r
 272   gr->GetXaxis()->SetTitle("Tan(#theta)");\r
 273   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV/c)}");\r
 274   gr->GetYaxis()->SetTitle("#sigmaDCA (cm)");\r
 275   gr->GetHistogram()->Write("DCAResolSPTTan");\r
 276 \r
 277   gPad->Clear();\r
 278   gr0->Draw("al*");\r
 279 \r
 280   fp->Close();\r
 281 }\r