PWG1/AliComparisonDCA.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliComparisonDCA class. It keeps information from
   3 // comparison of reconstructed and MC particle tracks. In addtion,
   4 // it keeps selection cuts used during comparison. The comparison
   5 // information is stored in the ROOT histograms. Analysis of these
   6 // histograms can be done by using Analyse() class function. The result of
   7 // the analysis (histograms/graphs) are stored in the folder
   8 // which is a data member of AliComparisonDCA.
   9 //
  10 // Author: J.Otwinowski 04/02/2008
  11 //------------------------------------------------------------------------------
  12
  13 /*
  14
  15   // after running comparison task, read the file, and get component
  16   gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
  17   LoadMyLibs();
  18   TFile f("Output.root");
  19   AliComparisonDCA * compObj = (AliComparisonDCA*)f.Get("AliComparisonDCA");
  20
  21   // Analyse comparison data
  22   compObj->Analyse();
  23
  24   // the output histograms/graphs will be stored in the folder "folderDCA"
  25   compObj->GetAnalysisFolder()->ls("*");
  26
  27   // user can save whole comparison object (or only folder with anlysed histograms)
  28   // in the seperate output file (e.g.)
  29   TFile fout("Analysed_DCA.root","recreate");
  30   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  31   fout.Close();
  32
  33 */
  34
  35 #include <iostream>
  36
  37 #include "TFile.h"
  38 #include "TCint.h"
  39 #include "TH3F.h"
  40 #include "TH2F.h"
  41 #include "TF1.h"
  42 #include "TProfile.h"
  43 #include "TProfile2D.h"
  44 #include "TGraph2D.h"
  45 #include "TCanvas.h"
  46 #include "TGraph.h"
  47 //
  48 #include "AliTracker.h"
  49 #include "AliESDEvent.h"
  50 #include "AliESD.h"
  51 #include "AliESDfriend.h"
  52 #include "AliESDfriendTrack.h"
  53 #include "AliRecInfoCuts.h"
  54 #include "AliMCInfoCuts.h"
  55 #include "AliLog.h"
  56 #include "AliESDVertex.h"
  57 //
  58 #include "AliMathBase.h"
  59 #include "AliTreeDraw.h"
  60
  61 #include "AliMCInfo.h"
  62 #include "AliESDRecInfo.h"
  63 #include "AliComparisonDCA.h"
  64
  65 using namespace std;
  66
  67 ClassImp(AliComparisonDCA)
  68
  69 //_____________________________________________________________________________
  70 AliComparisonDCA::AliComparisonDCA():
  71 //  TNamed("AliComparisonDCA","AliComparisonDCA"),
  72   AliComparisonObject("AliComparisonDCA"),
  73
  74   // DCA histograms
  75   fD0TanSPtB1(0),
  76   fD1TanSPtB1(0),
  77   fD0TanSPtL1(0),
  78   fD1TanSPtL1(0),
  79   fD0TanSPtInTPC(0),
  80   fD1TanSPtInTPC(0),
  81   fVertex(0),
  82
  83   // Cuts
  84   fCutsRC(0),
  85   fCutsMC(0),
  86
  87   // histogram folder
  88   fAnalysisFolder(0)
  89 {
  90   Init();
  91 }
  92
  93 //_____________________________________________________________________________
  94 AliComparisonDCA::~AliComparisonDCA()
  95 {
  96   //
  97   if(fD0TanSPtB1) delete fD0TanSPtB1; fD0TanSPtB1=0;
  98   if(fD1TanSPtB1) delete fD1TanSPtB1; fD1TanSPtB1=0;
  99   if(fD0TanSPtL1) delete fD0TanSPtL1; fD0TanSPtL1=0;
 100   if(fD1TanSPtL1) delete fD1TanSPtL1; fD1TanSPtL1=0;
 101   if(fD0TanSPtInTPC) delete fD0TanSPtInTPC; fD0TanSPtInTPC=0;
 102   if(fD1TanSPtInTPC) delete fD1TanSPtInTPC; fD1TanSPtInTPC=0;
 103   if(fVertex) delete fVertex; fVertex=0;
 104   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 105
 106 }
 107
 108 //_____________________________________________________________________________
 109 void AliComparisonDCA::Init()
 110 {
 111   // DCA histograms
 112   fD0TanSPtB1 = new TH3F("DCAyTanSPtB1","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
 113   fD0TanSPtB1->SetXTitle("tan(#theta)");
 114   fD0TanSPtB1->SetYTitle("#sqrt{p_{t}(GeV/c)}");
 115   fD0TanSPtB1->SetZTitle("DCA_{xy}");
 116
 117   fD1TanSPtB1 = new TH3F("DCAzTanSPtB1","DCAzTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
 118   fD1TanSPtB1->SetXTitle("tan(#theta)");
 119   fD1TanSPtB1->SetYTitle("#sqrt(p_{t}(GeV/c))");
 120   fD1TanSPtB1->SetZTitle("DCA_{z}");
 121
 122   fD0TanSPtL1 = new TH3F("DCAyTanSPtL1","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
 123   fD0TanSPtL1->SetXTitle("tan(#theta)");
 124   fD0TanSPtL1->SetYTitle("#sqrt{p_{t}(GeV/c)}");
 125   fD0TanSPtL1->SetZTitle("DCA_{xy}");
 126
 127   fD1TanSPtL1 = new TH3F("DCAzTanSPtL1","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1);
 128   fD1TanSPtL1->SetXTitle("tan(#theta)");
 129   fD1TanSPtL1->SetYTitle("#sqrt{p_{t}(GeV/c)}");
 130   fD1TanSPtL1->SetZTitle("DCA_{z}");
 131
 132   fD0TanSPtInTPC = new TH3F("DCAyTanSPtInTPC","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
 133   fD0TanSPtInTPC->SetXTitle("tan(#theta)");
 134   fD0TanSPtInTPC->SetYTitle("#sqrt{p_{t}(GeV/c)}");
 135   fD0TanSPtInTPC->SetZTitle("DCA_{xy}");
 136
 137   fD1TanSPtInTPC = new TH3F("DCAzTanSPtInTPC","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1);
 138   fD1TanSPtInTPC->SetXTitle("tan(#theta)");
 139   fD1TanSPtInTPC->SetYTitle("#sqrt{p_{t}(GeV/c)}");
 140   fD1TanSPtInTPC->SetZTitle("DCA_{z}");
 141
 142   // init cuts
 143   if(!fCutsMC)
 144     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 145   if(!fCutsRC)
 146     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 147
 148   // init folder
 149   fAnalysisFolder = CreateFolder("folderDCA","Analysis DCA Folder");
 150
 151   // vertex (0,0,0)
 152   fVertex = new AliESDVertex();
 153   fVertex->SetXv(0.0);
 154   fVertex->SetYv(0.0);
 155   fVertex->SetZv(0.0);
 156 }
 157
 158 //_____________________________________________________________________________
 159 void AliComparisonDCA::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
 160 {
 161   // Fill DCA comparison information
 162   AliExternalTrackParam *track = 0;
 163   Double_t field      = AliTracker::GetBz(); // nominal Bz field [kG]
 164   Double_t kMaxD      = 123456.0; // max distance
 165
 166   Int_t clusterITS[200];
 167   Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 168   Float_t dca1[2], cov1[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 169
 170   Float_t mcpt = infoMC->GetParticle().Pt();
 171   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
 172   Float_t spt = TMath::Sqrt(mcpt);
 173
 174   // distance to Prim. vertex
 175   const Double_t* dv = infoMC->GetVDist();
 176
 177   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
 178
 179   // Check selection cuts
 180   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 181   if (!isPrim) return;
 182   if (infoRC->GetStatus(1)!=3) return;
 183   if (!infoRC->GetESDtrack()) return;
 184   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 185   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
 186
 187   // calculate and set prim. vertex
 188   fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
 189   fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
 190   fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
 191
 192   // calculate track parameters at vertex
 193   if (infoRC->GetESDtrack()->GetTPCInnerParam())
 194   {
 195     if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
 196     {
 197       Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
 198
 199       if(bDCAStatus) {
 200         fD0TanSPtInTPC->Fill(tantheta,spt,dca[0]);
 201         fD1TanSPtInTPC->Fill(tantheta,spt,dca[1]);
 202           }
 203
 204           delete track;
 205     }
 206   }
 207
 208  // ITS + TPC
 209  infoRC->GetESDtrack()->GetImpactParameters(dca1,cov1);
 210
 211  if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)==0){
 212     fD0TanSPtB1->Fill(tantheta,spt,dca1[0]);
 213     fD1TanSPtB1->Fill(tantheta,spt,dca1[1]);
 214   }
 215     fD0TanSPtL1->Fill(tantheta,spt,dca1[0]);
 216     fD1TanSPtL1->Fill(tantheta,spt,dca1[1]);
 217 }
 218
 219 //_____________________________________________________________________________
 220 Long64_t AliComparisonDCA::Merge(TCollection* list)
 221 {
 222   // Merge list of objects (needed by PROOF)
 223
 224   if (!list)
 225   return 0;
 226
 227   if (list->IsEmpty())
 228   return 1;
 229
 230   TIterator* iter = list->MakeIterator();
 231   TObject* obj = 0;
 232
 233   // collection of generated histograms
 234   Int_t count=0;
 235   while((obj = iter->Next()) != 0)
 236   {
 237     AliComparisonDCA* entry = dynamic_cast<AliComparisonDCA*>(obj);
 238     if (entry == 0) continue;
 239
 240
 241     fD0TanSPtB1->Add(entry->fD0TanSPtB1);
 242     fD1TanSPtB1->Add(entry->fD1TanSPtB1);
 243     fD0TanSPtL1->Add(entry->fD0TanSPtL1);
 244     fD1TanSPtL1->Add(entry->fD1TanSPtL1);
 245     fD0TanSPtInTPC->Add(entry->fD0TanSPtInTPC);
 246     fD1TanSPtInTPC->Add(entry->fD1TanSPtInTPC);
 247
 248     count++;
 249   }
 250
 251 return count;
 252 }
 253
 254 //_____________________________________________________________________________
 255 void AliComparisonDCA::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
 256   // Process comparison information
 257   Process(infoMC,infoRC);
 258 }
 259
 260 //_____________________________________________________________________________
 261 void AliComparisonDCA::Analyse()
 262 {
 263   //
 264   // Analyse comparison information and store output histograms
 265   // in the analysis folder "folderDCA"
 266   //
 267
 268   TH1::AddDirectory(kFALSE);
 269
 270   TGraph * gr[4]= { 0,0,0,0 };
 271   TGraph2D *gr2[4]= { 0,0,0,0};
 272   AliComparisonDCA * comp=this;
 273   TObjArray *aFolderObj = new TObjArray;
 274
 275   // write results in the folder
 276   // Canvas to draw analysed histograms
 277   TCanvas * c = new TCanvas("canDCA","DCA resolution");
 278   c->Divide(2,4);
 279   //
 280   // DCA resolution
 281   //
 282   c->cd(1);
 283   gr[0] = AliMathBase::MakeStat1D(comp->fD0TanSPtB1,2,5);
 284   gr[0]->GetXaxis()->SetTitle("Tan(#theta)");
 285   gr[0]->GetYaxis()->SetTitle("#sigmaDCA_xy (cm)");
 286   gr[0]->SetName("DCAXYResolTan");
 287   gr[0]->Draw("Al*");
 288
 289   aFolderObj->Add(gr[0]);
 290
 291   c->cd(2);
 292   gr[1] = AliMathBase::MakeStat1D(comp->fD1TanSPtB1,2,5);
 293   gr[1]->GetXaxis()->SetTitle("Tan(#theta)");
 294   gr[1]->GetYaxis()->SetTitle("#sigmaDCA_z (cm)");
 295   gr[1]->SetName("DCAZResolTan");
 296   gr[1]->Draw("Al*");
 297
 298   aFolderObj->Add(gr[1]);
 299
 300   //
 301   // DCA mean value
 302   //
 303   c->cd(3);
 304   gr[2] = AliMathBase::MakeStat1D(comp->fD0TanSPtB1,2,4);
 305   gr[2]->GetXaxis()->SetTitle("Tan(#theta)");
 306   gr[2]->GetYaxis()->SetTitle("mean DCA_xy (cm)");
 307   gr[2]->SetName("DCAXYMeanTan");
 308   gr[2]->Draw("Al*");
 309
 310   aFolderObj->Add(gr[2]);
 311
 312   c->cd(4);
 313   gr[3] = AliMathBase::MakeStat1D(comp->fD1TanSPtB1,2,4);
 314   gr[3]->GetXaxis()->SetTitle("Tan(#theta)");
 315   gr[3]->GetYaxis()->SetTitle("mean DCA_z (cm)");
 316   gr[3]->SetName("DCAZMeanTan");
 317   gr[3]->Draw("Al*");
 318
 319   aFolderObj->Add(gr[3]);
 320
 321   // 2D DCA resolution
 322   c->cd(5);
 323   gr2[0] = AliMathBase::MakeStat2D(comp->fD0TanSPtB1,4,2,5);
 324   gr2[0]->GetXaxis()->SetTitle("Tan(#theta)");
 325   gr2[0]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV/c)}");
 326   gr2[0]->GetZaxis()->SetTitle("#sigmaDCA_xy (cm)");
 327   gr2[0]->SetName("DCAXYResolSPTTan");
 328   gr2[0]->GetHistogram()->Draw("colz");
 329
 330   gr2[0]->GetHistogram()->SetName("DCAXYResolSPTTan");
 331   aFolderObj->Add(gr2[0]->GetHistogram());
 332
 333   c->cd(6);
 334   gr2[1] = AliMathBase::MakeStat2D(comp->fD1TanSPtB1,4,2,5);
 335   gr2[1]->GetXaxis()->SetTitle("Tan(#theta)");
 336   gr2[1]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV/c)}");
 337   gr2[1]->GetZaxis()->SetTitle("#sigmaDCA_z (cm)");
 338   gr2[1]->SetName("DCAZResolSPTTan");
 339   gr2[1]->GetHistogram()->Draw("colz");
 340
 341   gr2[1]->GetHistogram()->SetName("DCAZResolSPTTan");
 342   aFolderObj->Add(gr2[1]->GetHistogram());
 343
 344   // 2D DCA mean value
 345   c->cd(7);
 346   gr2[2] = AliMathBase::MakeStat2D(comp->fD0TanSPtB1,4,2,4);
 347   gr2[2]->GetXaxis()->SetTitle("Tan(#theta)");
 348   gr2[2]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV/c)}");
 349   gr2[2]->GetZaxis()->SetTitle("mean DCA_xy (cm)");
 350   gr2[2]->SetName("DCAXYMeanSPTTan");
 351   gr2[2]->GetHistogram()->Draw("colz");
 352
 353   gr2[2]->GetHistogram()->SetName("DCAXYMeanSPTTan");
 354   aFolderObj->Add(gr2[2]->GetHistogram());
 355
 356   c->cd(8);
 357   gr2[3] = AliMathBase::MakeStat2D(comp->fD1TanSPtB1,4,2,4);
 358   gr2[3]->GetXaxis()->SetTitle("Tan(#theta)");
 359   gr2[3]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV/c)}");
 360   gr2[3]->GetZaxis()->SetTitle("mean DCA_z (cm)");
 361   gr2[3]->SetName("DCAZMeanSPTTan");
 362   gr2[3]->GetHistogram()->Draw("colz");
 363
 364   gr2[3]->GetHistogram()->SetName("DCAZMeanSPTTan");
 365   aFolderObj->Add(gr2[3]->GetHistogram());
 366
 367   // export objects to analysis folder
 368   fAnalysisFolder = ExportToFolder(aFolderObj);
 369
 370   // delete only TObjArray
 371   if(aFolderObj) delete aFolderObj;
 372 }
 373
 374 //_____________________________________________________________________________
 375 TFolder* AliComparisonDCA::ExportToFolder(TObjArray * array)
 376 {
 377   // recreate folder avery time and export objects to new one
 378   //
 379   AliComparisonDCA * comp=this;
 380   TFolder *folder = comp->GetAnalysisFolder();
 381
 382   TString name, title;
 383   TFolder *newFolder = 0;
 384   Int_t i = 0;
 385   Int_t size = array->GetSize();
 386
 387   if(folder) {
 388      // get name and title from old folder
 389      name = folder->GetName();
 390      title = folder->GetTitle();
 391
 392          // delete old one
 393      delete folder;
 394
 395          // create new one
 396      newFolder = CreateFolder(name.Data(),title.Data());
 397      newFolder->SetOwner();
 398
 399          // add objects to folder
 400      while(i < size) {
 401            newFolder->Add(array->At(i));
 402            i++;
 403          }
 404   }
 405
 406 return newFolder;
 407 }
 408
 409
 410 //_____________________________________________________________________________
 411 TFolder* AliComparisonDCA::CreateFolder(TString name,TString title) {
 412 // create folder for analysed histograms
 413 TFolder *folder = 0;
 414   folder = new TFolder(name.Data(),title.Data());
 415
 416   return folder;
 417 }