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
 169   Float_t mcpt = infoMC->GetParticle().Pt();
 170   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
 171   Float_t spt = TMath::Sqrt(mcpt);
 172
 173   // distance to Prim. vertex
 174   const Double_t* dv = infoMC->GetVDist();
 175
 176   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
 177
 178   // Check selection cuts
 179   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 180   if (!isPrim) return;
 181   if (infoRC->GetStatus(1)!=3) return;
 182   if (!infoRC->GetESDtrack()) return;
 183   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 184   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
 185
 186   // calculate and set prim. vertex
 187   fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
 188   fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
 189   fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
 190
 191   // calculate track parameters at vertex
 192   if (infoRC->GetESDtrack()->GetTPCInnerParam())
 193   {
 194     if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
 195     {
 196       Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
 197
 198       if(bDCAStatus) {
 199         fD0TanSPtInTPC->Fill(tantheta,spt,dca[0]);
 200         fD1TanSPtInTPC->Fill(tantheta,spt,dca[1]);
 201           }
 202
 203           delete track;
 204     }
 205   }
 206
 207  if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)==0){
 208     fD0TanSPtB1->Fill(tantheta,spt,dca[0]);
 209     fD1TanSPtB1->Fill(tantheta,spt,dca[1]);
 210   }
 211     fD0TanSPtL1->Fill(tantheta,spt,dca[0]);
 212     fD1TanSPtL1->Fill(tantheta,spt,dca[1]);
 213 }
 214
 215 //_____________________________________________________________________________
 216 Long64_t AliComparisonDCA::Merge(TCollection* list)
 217 {
 218   // Merge list of objects (needed by PROOF)
 219
 220   if (!list)
 221   return 0;
 222
 223   if (list->IsEmpty())
 224   return 1;
 225
 226   TIterator* iter = list->MakeIterator();
 227   TObject* obj = 0;
 228
 229   // collection of generated histograms
 230   Int_t count=0;
 231   while((obj = iter->Next()) != 0)
 232   {
 233     AliComparisonDCA* entry = dynamic_cast<AliComparisonDCA*>(obj);
 234     if (entry == 0) continue;
 235
 236
 237     fD0TanSPtB1->Add(entry->fD0TanSPtB1);
 238     fD1TanSPtB1->Add(entry->fD1TanSPtB1);
 239     fD0TanSPtL1->Add(entry->fD0TanSPtL1);
 240     fD1TanSPtL1->Add(entry->fD1TanSPtL1);
 241     fD0TanSPtInTPC->Add(entry->fD0TanSPtInTPC);
 242     fD1TanSPtInTPC->Add(entry->fD1TanSPtInTPC);
 243
 244     count++;
 245   }
 246
 247 return count;
 248 }
 249
 250 //_____________________________________________________________________________
 251 void AliComparisonDCA::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
 252   // Process comparison information
 253   Process(infoMC,infoRC);
 254 }
 255
 256 //_____________________________________________________________________________
 257 void AliComparisonDCA::Analyse()
 258 {
 259   //
 260   // Analyse comparison information and store output histograms
 261   // in the analysis folder "folderDCA"
 262   //
 263
 264   TH1::AddDirectory(kFALSE);
 265
 266   TGraph2D *gr=0;
 267   TGraph * gr0=0;
 268   AliComparisonDCA * comp=this;
 269   TObjArray *aFolderObj = new TObjArray;
 270
 271   // write results in the folder
 272   // Canvas to draw analysed histograms
 273   TCanvas * c = new TCanvas("canDCA","DCA resolution");
 274   c->Divide(1,2);
 275   c->cd(1);
 276   //
 277   // DCA resolution
 278   //
 279   gr0 = AliMathBase::MakeStat1D(comp->fD0TanSPtB1,2,5);
 280   gr0->GetXaxis()->SetTitle("Tan(#theta)");
 281   gr0->GetYaxis()->SetTitle("#sigmaDCA (cm)");
 282   gr0->SetName("DCAResolTan");
 283   gr0->Draw("Al*");
 284
 285   //if(folder) folder->Add(gr0->GetHistogram());
 286   aFolderObj->Add(gr0);
 287   //
 288   c->cd(2);
 289   gr = AliMathBase::MakeStat2D(comp->fD0TanSPtB1,4,2,5);
 290   gr->GetXaxis()->SetTitle("Tan(#theta)");
 291   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV/c)}");
 292   gr->GetYaxis()->SetTitle("#sigmaDCA (cm)");
 293   gr->SetName("DCAResolSPTTan");
 294   gr->GetHistogram()->Draw("colz");
 295
 296   aFolderObj->Add(gr->GetHistogram());
 297
 298   // export objects to analysis folder
 299   fAnalysisFolder = ExportToFolder(aFolderObj);
 300
 301   // delete only TObjArray
 302   if(aFolderObj) delete aFolderObj;
 303 }
 304
 305 //_____________________________________________________________________________
 306 TFolder* AliComparisonDCA::ExportToFolder(TObjArray * array)
 307 {
 308   // recreate folder avery time and export objects to new one
 309   //
 310   AliComparisonDCA * comp=this;
 311   TFolder *folder = comp->GetAnalysisFolder();
 312
 313   TString name, title;
 314   TFolder *newFolder = 0;
 315   Int_t i = 0;
 316   Int_t size = array->GetSize();
 317
 318   if(folder) {
 319      // get name and title from old folder
 320      name = folder->GetName();
 321      title = folder->GetTitle();
 322
 323          // delete old one
 324      delete folder;
 325
 326          // create new one
 327      newFolder = CreateFolder(name.Data(),title.Data());
 328      newFolder->SetOwner();
 329
 330          // add objects to folder
 331      while(i < size) {
 332            newFolder->Add(array->At(i));
 333            i++;
 334          }
 335   }
 336
 337 return newFolder;
 338 }
 339
 340
 341 //_____________________________________________________________________________
 342 TFolder* AliComparisonDCA::CreateFolder(TString name,TString title) {
 343 // create folder for analysed histograms
 344 TFolder *folder = 0;
 345   folder = new TFolder(name.Data(),title.Data());
 346
 347   return folder;
 348 }