PWG1/AliComparisonDEdx.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliComparisonDEdx 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 which is
   8 // a data of AliComparisonDEdx.
   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   AliComparisonDEdx * compObj = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");
  20
  21   // Analyse comparison data
  22   compObj->Analyse();
  23
  24   // the output histograms/graphs will be stored in the folder "folderDEdx"
  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_DEdx.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 "AliESDEvent.h"
  49 #include "AliESD.h"
  50 #include "AliESDfriend.h"
  51 #include "AliESDfriendTrack.h"
  52 #include "AliRecInfoCuts.h"
  53 #include "AliMCInfoCuts.h"
  54 #include "AliLog.h"
  55 //
  56 #include "AliMathBase.h"
  57 #include "AliTreeDraw.h"
  58 //#include "TStatToolkit.h"
  59
  60 #include "AliMCInfo.h"
  61 #include "AliESDRecInfo.h"
  62 #include "AliComparisonDEdx.h"
  63
  64 using namespace std;
  65
  66 ClassImp(AliComparisonDEdx)
  67
  68 //_____________________________________________________________________________
  69 AliComparisonDEdx::AliComparisonDEdx():
  70 //  TNamed("AliComparisonDEdx","AliComparisonDEdx"),
  71   AliComparisonObject("AliComparisonDEdx"),
  72
  73   // dEdx
  74   fTPCSignalNormTan(0),
  75   fTPCSignalNormSPhi(0),
  76   fTPCSignalNormTPhi(0),
  77   //
  78   fTPCSignalNormTanSPhi(0),
  79   fTPCSignalNormTanTPhi(0),
  80   fTPCSignalNormTanSPt(0),
  81
  82   // Cuts
  83   fCutsRC(0),
  84   fCutsMC(0),
  85   fMCPtMin(0),
  86   fMCAbsTanThetaMax(0),
  87   fMCPdgCode(0),
  88
  89   // histogram folder
  90   fAnalysisFolder(0)
  91 {
  92   Init();
  93 }
  94
  95 //_____________________________________________________________________________
  96 AliComparisonDEdx::~AliComparisonDEdx(){
  97
  98   if(fTPCSignalNormTan)  delete fTPCSignalNormTan; fTPCSignalNormTan=0;
  99   if(fTPCSignalNormSPhi) delete fTPCSignalNormSPhi; fTPCSignalNormSPhi=0;
 100   if(fTPCSignalNormTPhi) delete fTPCSignalNormTPhi; fTPCSignalNormTPhi=0;
 101   //
 102   if(fTPCSignalNormTanSPhi) delete fTPCSignalNormTanSPhi; fTPCSignalNormTanSPhi=0;
 103   if(fTPCSignalNormTanTPhi) delete fTPCSignalNormTanTPhi; fTPCSignalNormTanTPhi=0;
 104   if(fTPCSignalNormTanSPt)  delete fTPCSignalNormTanSPt; fTPCSignalNormTanSPt=0;
 105
 106   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 107 }
 108
 109 //_____________________________________________________________________________
 110 void AliComparisonDEdx::Init()
 111 {
 112   // Init histograms
 113
 114   // TPC dEdx
 115   fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2,  40,30,70);
 116   fTPCSignalNormTan->SetXTitle("tan(#theta)");
 117   fTPCSignalNormTan->SetYTitle("rec. dE/dx / calc. dE/dx");
 118
 119   fTPCSignalNormSPhi   = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70);
 120   fTPCSignalNormSPhi->SetXTitle("sin(#phi)");
 121   fTPCSignalNormSPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
 122
 123   fTPCSignalNormTPhi   = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70);
 124   fTPCSignalNormTPhi->SetXTitle("tan(#phi)");
 125   fTPCSignalNormTPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
 126
 127   fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
 128   fTPCSignalNormTanSPhi->SetXTitle("tan(#theta)");
 129   fTPCSignalNormTanSPhi->SetYTitle("sin(#phi)");
 130   fTPCSignalNormTanSPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
 131
 132   fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
 133   fTPCSignalNormTanTPhi->SetXTitle("tan(#theta)");
 134   fTPCSignalNormTanTPhi->SetYTitle("tan(#phi)");
 135   fTPCSignalNormTanTPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
 136
 137   fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3, 40,30,70);
 138   fTPCSignalNormTanSPt->SetXTitle("tan(#theta)");
 139   fTPCSignalNormTanSPt->SetYTitle("#sqrt{p_{t}}");
 140   fTPCSignalNormTanSPt->SetZTitle("rec. dE/dx / calc. dE/dx");
 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("folderDEdx","Analysis de/dx Folder");
 150 }
 151
 152 //_____________________________________________________________________________
 153 void AliComparisonDEdx::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
 154
 155   // Fill dE/dx  comparison information
 156
 157   Float_t mcpt = infoMC->GetParticle().Pt();
 158   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
 159   Float_t mprim = infoMC->GetPrim();
 160
 161   // distance to Prim. vertex
 162   const Double_t* dv = infoMC->GetVDist();
 163
 164   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
 165
 166   // Check selection cuts
 167   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 168   if (!isPrim) return;
 169   if (infoRC->GetStatus(1)!=3) return;
 170   if (!infoRC->GetESDtrack()) return;
 171   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 172   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
 173   //if (mprim>1.4) return;
 174   //if (mprim<0.5) return;
 175   if (mprim > fCutsMC->GetMaxTPCSignal()) return;
 176   if (mprim < fCutsMC->GetMinTPCSignal()) return;
 177   if (infoRC->GetESDtrack()->GetTPCsignalN()<fCutsRC->GetMinTPCsignalN()) return;
 178   //
 179   Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();
 180   Float_t sphi =  infoRC->GetESDtrack()->GetInnerParam()->GetSnp();
 181   Float_t tphi =  sphi/TMath::Sqrt(1-sphi*sphi);
 182
 183   if (TMath::Abs(infoMC->GetParticle().GetPdgCode()) != GetMCPdgCode()) return;
 184   //if (mcpt>0.5) {
 185   if (mcpt > GetMCPtMin()) {
 186     fTPCSignalNormTan->Fill(tantheta,ratio);    // only subset
 187   }
 188
 189   //if (TMath::Abs(tantheta)<0.5){
 190   if (TMath::Abs(tantheta) < GetMCAbsTanThetaMax()){
 191     fTPCSignalNormSPhi->Fill(sphi,ratio);       // only subset
 192     fTPCSignalNormTPhi->Fill(tphi,ratio);       // only subset
 193   }
 194   fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio);
 195   fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio);
 196   fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);
 197 }
 198
 199 //_____________________________________________________________________________
 200 Long64_t AliComparisonDEdx::Merge(TCollection* list)
 201 {
 202   // Merge list of objects (needed by PROOF)
 203
 204   if (!list)
 205   return 0;
 206
 207   if (list->IsEmpty())
 208   return 1;
 209
 210   TIterator* iter = list->MakeIterator();
 211   TObject* obj = 0;
 212
 213   // collection of generated histograms
 214   Int_t count=0;
 215   while((obj = iter->Next()) != 0)
 216   {
 217     AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);
 218     if (entry == 0) continue;
 219
 220     fTPCSignalNormTan->Add(entry->fTPCSignalNormTan);
 221     fTPCSignalNormSPhi->Add(entry->fTPCSignalNormSPhi);
 222     fTPCSignalNormTPhi->Add(entry->fTPCSignalNormTPhi);
 223     //
 224     fTPCSignalNormTanSPhi->Add(entry->fTPCSignalNormTanSPhi);
 225     fTPCSignalNormTanTPhi->Add(entry->fTPCSignalNormTanTPhi);
 226     fTPCSignalNormTanSPt->Add(entry->fTPCSignalNormTanSPt);
 227
 228     count++;
 229   }
 230
 231 return count;
 232 }
 233
 234 //_____________________________________________________________________________
 235 void AliComparisonDEdx::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
 236 {
 237   // Process comparison information
 238   Process(infoMC,infoRC);
 239 }
 240
 241 //_____________________________________________________________________________
 242 TH1F* AliComparisonDEdx::MakeResol(TH2F * his, Int_t integ, Bool_t type)
 243 {
 244   // Make resolution histograms
 245   TH1F *hisr, *hism;
 246   if (!gPad) new TCanvas;
 247   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
 248   if (type) return hism;
 249   else
 250     return hisr;
 251 }
 252
 253 //_____________________________________________________________________________
 254 void AliComparisonDEdx::Analyse()
 255 {
 256   // Analyze comparison information and store output histograms
 257   // in the folder "folderDEdx"
 258   //
 259
 260   TH1::AddDirectory(kFALSE);
 261
 262   AliComparisonDEdx * comp=this;
 263   TObjArray *aFolderObj = new TObjArray;
 264
 265   TH1F *hiss=0;
 266   TGraph2D * gr=0;
 267
 268   // write results in the folder
 269   TCanvas * c = new TCanvas("can","TPC dedx");
 270   c->cd();
 271
 272   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);
 273   hiss->SetXTitle("Tan(#theta)");
 274   hiss->SetYTitle("#sigma_{dEdx}");
 275   hiss->Draw();
 276   hiss->SetName("TPCdEdxResolTan");
 277
 278   aFolderObj->Add(hiss);
 279   //
 280   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1);
 281   hiss->SetXTitle("Tan(#theta)");
 282   hiss->SetYTitle("<dEdx>");
 283   hiss->Draw();
 284   hiss->SetName("TPCdEdxMeanTan");
 285
 286   aFolderObj->Add(hiss);
 287   //
 288   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);
 289   gr->GetXaxis()->SetTitle("Tan(#theta)");
 290   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
 291   gr->GetZaxis()->SetTitle("<dEdx>");
 292   gr->SetName("TPCdEdxMeanTanPt_1");
 293   gr->GetHistogram()->Draw("colz");
 294
 295   gr->GetHistogram()->SetName("TPCdEdxMeanTanPt_1");
 296   aFolderObj->Add(gr->GetHistogram());
 297   //
 298   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);
 299   gr->GetXaxis()->SetTitle("Tan(#theta)");
 300   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
 301   gr->GetZaxis()->SetTitle("#sigma_{dEdx}");
 302   gr->SetName("TPCdEdxMeanTanPt_2");
 303   gr->GetHistogram()->Draw("colz");
 304
 305   gr->GetHistogram()->SetName("TPCdEdxMeanTanPt_2");
 306   aFolderObj->Add(gr->GetHistogram());
 307
 308   // export objects to analysis folder
 309   fAnalysisFolder = ExportToFolder(aFolderObj);
 310
 311   // delete only TObjrArray
 312   if(aFolderObj) delete aFolderObj;
 313 }
 314
 315 //_____________________________________________________________________________
 316 TFolder* AliComparisonDEdx::ExportToFolder(TObjArray * array)
 317 {
 318   // recreate folder avery time and export objects to new one
 319   //
 320   AliComparisonDEdx * comp=this;
 321   TFolder *folder = comp->GetAnalysisFolder();
 322
 323   TString name, title;
 324   TFolder *newFolder = 0;
 325   Int_t i = 0;
 326   Int_t size = array->GetSize();
 327
 328   if(folder) {
 329      // get name and title from old folder
 330      name = folder->GetName();
 331      title = folder->GetTitle();
 332
 333          // delete old one
 334      delete folder;
 335
 336          // create new one
 337      newFolder = CreateFolder(name.Data(),title.Data());
 338      newFolder->SetOwner();
 339
 340          // add objects to folder
 341      while(i < size) {
 342            newFolder->Add(array->At(i));
 343            i++;
 344          }
 345   }
 346
 347 return newFolder;
 348 }
 349
 350
 351 //_____________________________________________________________________________
 352 TFolder* AliComparisonDEdx::CreateFolder(TString name,TString title) {
 353 // create folder for analysed histograms
 354 TFolder *folder = 0;
 355   folder = new TFolder(name.Data(),title.Data());
 356
 357   return folder;
 358 }
 359