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   AliComparisonDEdx * compObj = (AliComparisonDEdx*)cOutput->FindObject("AliComparisonDEdx");
  21
  22   // Analyse comparison data
  23   compObj->Analyse();
  24
  25   // the output histograms/graphs will be stored in the folder "folderDEdx"
  26   compObj->GetAnalysisFolder()->ls("*");
  27
  28   // user can save whole comparison object (or only folder with anlysed histograms)
  29   // in the seperate output file (e.g.)
  30   TFile fout("Analysed_DEdx.root"."recreate");
  31   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  32   fout.Close();
  33
  34 */
  35
  36 #include <TAxis.h>
  37 #include <TCanvas.h>
  38 #include <TH1.h>
  39
  40 #include "AliComparisonDEdx.h"
  41 #include "AliESDEvent.h"
  42 #include "AliESDRecInfo.h"
  43 #include "AliLog.h"
  44 #include "AliMCInfo.h"
  45 #include "AliMCInfoCuts.h"
  46 #include "AliMathBase.h"
  47 #include "AliRecInfoCuts.h"
  48 #include "AliTreeDraw.h"
  49
  50 using namespace std;
  51
  52 ClassImp(AliComparisonDEdx)
  53
  54 //_____________________________________________________________________________
  55 AliComparisonDEdx::AliComparisonDEdx():
  56 //  TNamed("AliComparisonDEdx","AliComparisonDEdx"),
  57   AliComparisonObject("AliComparisonDEdx"),
  58
  59   // dEdx
  60   fDeDxHisto(0),
  61
  62   // Cuts
  63   fCutsRC(0),
  64   fCutsMC(0),
  65
  66   // histogram folder
  67   fAnalysisFolder(0)
  68 {
  69   // default constructor
  70 }
  71
  72 //_____________________________________________________________________________
  73 AliComparisonDEdx::AliComparisonDEdx(Char_t* name="AliComparisonDEdx", Char_t* title="AliComparisonDEdx",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
  74   AliComparisonObject(name,title),
  75
  76   // dEdx
  77   fDeDxHisto(0),
  78
  79   // Cuts
  80   fCutsRC(0),
  81   fCutsMC(0),
  82
  83   // histogram folder
  84   fAnalysisFolder(0)
  85 {
  86   // named constructor
  87
  88   SetAnalysisMode(analysisMode);
  89   SetHptGenerator(hptGenerator);
  90   Init();
  91 }
  92
  93
  94 //_____________________________________________________________________________
  95 AliComparisonDEdx::~AliComparisonDEdx()
  96 {
  97   // destructor
  98   if(fDeDxHisto)  delete fDeDxHisto; fDeDxHisto=0;
  99
 100   if(fCutsRC) delete fCutsRC; fCutsRC=0;
 101   if(fCutsMC) delete fCutsMC; fCutsMC=0;
 102
 103   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 104 }
 105
 106 //_____________________________________________________________________________
 107 void AliComparisonDEdx::Init()
 108 {
 109   // Init histograms
 110
 111   // TPC dEdx
 112   Int_t nPBins = 31;
 113     Double_t binsP[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
 114     Double_t pMin = 0., pMax = 10.;
 115
 116     if(IsHptGenerator() == kTRUE) {
 117       nPBins = 100;
 118       pMin = 0.; pMax = 100.;
 119     }
 120
 121    //signal:alpha:y:z:snp:tgl:ncls:pid:p
 122    Int_t binsQA[9]    = {600,50, 50,  50, 50, 50, 80, 5, nPBins};
 123    Double_t xminQA[9] = {0,  -4,-20,-250, -1, -2, 0, 0., pMin};
 124    Double_t xmaxQA[9] = {300, 4, 20, 250,  1,  2, 160, 5., pMax};
 125
 126    fDeDxHisto = new THnSparseF("fDeDxHisto","signal:alpha:y:z:snp:tgl:ncls:pid:momentum",9,binsQA,xminQA,xmaxQA);
 127    if(!IsHptGenerator()) fDeDxHisto->SetBinEdges(8,binsP);
 128
 129    fDeDxHisto->GetAxis(0)->SetTitle("signal");
 130    fDeDxHisto->GetAxis(1)->SetTitle("alpha (rad)");
 131    fDeDxHisto->GetAxis(2)->SetTitle("y (cm)");
 132    fDeDxHisto->GetAxis(3)->SetTitle("z (cm)");
 133    fDeDxHisto->GetAxis(4)->SetTitle("snp");
 134    fDeDxHisto->GetAxis(5)->SetTitle("tgl");
 135    fDeDxHisto->GetAxis(6)->SetTitle("ncls");
 136    fDeDxHisto->GetAxis(6)->SetTitle("pid");
 137    fDeDxHisto->GetAxis(7)->SetTitle("p (GeV/c)");
 138    fDeDxHisto->Sumw2();
 139
 140    // Init cuts
 141    if(!fCutsMC)
 142      AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 143    if(!fCutsRC)
 144      AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 145
 146    // init folder
 147    fAnalysisFolder = CreateFolder("folderDEdx","Analysis de/dx Folder");
 148 }
 149
 150 //_____________________________________________________________________________
 151 void AliComparisonDEdx::ProcessTPC(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
 152 {
 153   // Fill dE/dx  comparison information
 154
 155   // Check selection cuts
 156   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 157
 158   Double_t pid = -1;
 159   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetEM() ) pid = 0;
 160   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetMuM() ) pid = 1;
 161   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetPiP() ) pid = 2;
 162   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetKP() ) pid = 3;
 163   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetProt() ) pid = 4;
 164
 165   if (!infoRC->GetESDtrack()) return;
 166   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 167
 168   Float_t dedx = infoRC->GetESDtrack()->GetTPCsignal();
 169   Int_t ncls = infoRC->GetESDtrack()->GetTPCNcls();
 170
 171   const AliExternalTrackParam *innerParam =  0;
 172   if ((innerParam = infoRC->GetESDtrack()->GetInnerParam()) == 0) return;
 173
 174   Double_t pt = innerParam->Pt();
 175   Double_t lam = TMath::ATan2(innerParam->Pz(),innerParam->Pt());
 176   Double_t p = pt/TMath::Cos(lam);
 177   Double_t alpha = innerParam->GetAlpha();
 178   Double_t y = innerParam->GetY();
 179   Double_t z = innerParam->GetZ();
 180   Double_t snp = innerParam->GetSnp();
 181   Double_t tgl = innerParam->GetTgl();
 182
 183   Double_t vDeDxHisto[9] = {dedx,alpha,y,z,snp,tgl,ncls,pid,p};
 184   fDeDxHisto->Fill(vDeDxHisto);
 185 }
 186
 187 //_____________________________________________________________________________
 188 void AliComparisonDEdx::ProcessTPCITS(AliMCInfo* const /*infoMC*/, AliESDRecInfo *const /*infoRC*/)
 189 {
 190   // Fill dE/dx  comparison information
 191
 192    AliDebug(AliLog::kWarning, "Warning: Not implemented");
 193 }
 194
 195 //_____________________________________________________________________________
 196 void AliComparisonDEdx::ProcessConstrained(AliMCInfo* const /*infoMC*/, AliESDRecInfo *const /*infoRC*/)
 197 {
 198   // Fill dE/dx  comparison information
 199
 200    AliDebug(AliLog::kWarning, "Warning: Not implemented");
 201 }
 202
 203 //_____________________________________________________________________________
 204 Long64_t AliComparisonDEdx::Merge(TCollection* const list)
 205 {
 206   // Merge list of objects (needed by PROOF)
 207
 208   if (!list)
 209   return 0;
 210
 211   if (list->IsEmpty())
 212   return 1;
 213
 214   TIterator* iter = list->MakeIterator();
 215   TObject* obj = 0;
 216
 217   // collection of generated histograms
 218   Int_t count=0;
 219   while((obj = iter->Next()) != 0)
 220   {
 221     AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);
 222     if (entry == 0) continue;
 223
 224     fDeDxHisto->Add(entry->fDeDxHisto);
 225     count++;
 226   }
 227
 228 return count;
 229 }
 230
 231 //_____________________________________________________________________________
 232 void AliComparisonDEdx::Exec(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
 233 {
 234   // Process comparison information
 235
 236   if(GetAnalysisMode() == 0) ProcessTPC(infoMC,infoRC);
 237   else if(GetAnalysisMode() == 1) ProcessTPCITS(infoMC,infoRC);
 238   else if(GetAnalysisMode() == 2) ProcessConstrained(infoMC,infoRC);
 239   else {
 240     printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 241     return;
 242   }
 243 }
 244
 245 //_____________________________________________________________________________
 246 TH1F* AliComparisonDEdx::MakeResol(TH2F * his, Int_t integ, Bool_t type)
 247 {
 248   // Make resolution histograms
 249   TH1F *hisr, *hism;
 250   if (!gPad) new TCanvas;
 251   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
 252   if (type) return hism;
 253   else
 254     return hisr;
 255 }
 256
 257 //_____________________________________________________________________________
 258 void AliComparisonDEdx::Analyse()
 259 {
 260   // Analyze comparison information and store output histograms
 261   // in the folder "folderDEdx"
 262   //
 263   TH1::AddDirectory(kFALSE);
 264   TObjArray *aFolderObj = new TObjArray;
 265
 266   /*
 267   TH1F *hiss=0;
 268   TGraph2D * gr=0;
 269
 270   // write results in the folder
 271   TCanvas * c = new TCanvas("can","TPC dedx");
 272   c->cd();
 273
 274   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);
 275   hiss->SetXTitle("Tan(#theta)");
 276   hiss->SetYTitle("#sigma_{dEdx}");
 277   hiss->Draw();
 278   hiss->SetName("TPCdEdxResolTan");
 279
 280   aFolderObj->Add(hiss);
 281   //
 282   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1);
 283   hiss->SetXTitle("Tan(#theta)");
 284   hiss->SetYTitle("<dEdx>");
 285   hiss->Draw();
 286   hiss->SetName("TPCdEdxMeanTan");
 287
 288   aFolderObj->Add(hiss);
 289   //
 290   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);
 291   gr->GetXaxis()->SetTitle("Tan(#theta)");
 292   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
 293   gr->GetZaxis()->SetTitle("<dEdx>");
 294   gr->SetName("TPCdEdxMeanTanPt_1");
 295   gr->GetHistogram()->Draw("colz");
 296
 297   gr->GetHistogram()->SetName("TPCdEdxMeanTanPt_1");
 298   aFolderObj->Add(gr->GetHistogram());
 299   //
 300   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);
 301   gr->GetXaxis()->SetTitle("Tan(#theta)");
 302   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
 303   gr->GetZaxis()->SetTitle("#sigma_{dEdx}");
 304   gr->SetName("TPCdEdxMeanTanPt_2");
 305   gr->GetHistogram()->Draw("colz");
 306
 307   gr->GetHistogram()->SetName("TPCdEdxMeanTanPt_2");
 308   aFolderObj->Add(gr->GetHistogram());
 309   */
 310
 311   // export objects to analysis folder
 312   fAnalysisFolder = ExportToFolder(aFolderObj);
 313
 314   // delete only TObjrArray
 315   if(aFolderObj) delete aFolderObj;
 316 }
 317
 318 //_____________________________________________________________________________
 319 TFolder* AliComparisonDEdx::ExportToFolder(TObjArray * array)
 320 {
 321   // recreate folder avery time and export objects to new one
 322   //
 323   AliComparisonDEdx * comp=this;
 324   TFolder *folder = comp->GetAnalysisFolder();
 325
 326   TString name, title;
 327   TFolder *newFolder = 0;
 328   Int_t i = 0;
 329   Int_t size = array->GetSize();
 330
 331   if(folder) {
 332      // get name and title from old folder
 333      name = folder->GetName();
 334      title = folder->GetTitle();
 335
 336          // delete old one
 337      delete folder;
 338
 339          // create new one
 340      newFolder = CreateFolder(name.Data(),title.Data());
 341      newFolder->SetOwner();
 342
 343          // add objects to folder
 344      while(i < size) {
 345            newFolder->Add(array->At(i));
 346            i++;
 347          }
 348   }
 349
 350 return newFolder;
 351 }
 352
 353
 354 //_____________________________________________________________________________
 355 TFolder* AliComparisonDEdx::CreateFolder(TString name,TString title) {
 356 // create folder for analysed histograms
 357 TFolder *folder = 0;
 358   folder = new TFolder(name.Data(),title.Data());
 359
 360   return folder;
 361 }
 362