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 "TCanvas.h"
  37 #include "AliESDEvent.h"
  38 #include "AliRecInfoCuts.h"
  39 #include "AliMCInfoCuts.h"
  40 #include "AliLog.h"
  41 #include "AliMathBase.h"
  42 #include "AliTreeDraw.h"
  43
  44 #include "AliMCInfo.h"
  45 #include "AliESDRecInfo.h"
  46 #include "AliComparisonDEdx.h"
  47
  48 using namespace std;
  49
  50 ClassImp(AliComparisonDEdx)
  51
  52 //_____________________________________________________________________________
  53 AliComparisonDEdx::AliComparisonDEdx():
  54 //  TNamed("AliComparisonDEdx","AliComparisonDEdx"),
  55   AliComparisonObject("AliComparisonDEdx"),
  56
  57   // dEdx
  58   fDeDxHisto(0),
  59
  60   // Cuts
  61   fCutsRC(0),
  62   fCutsMC(0),
  63
  64   // histogram folder
  65   fAnalysisFolder(0)
  66 {
  67   // default constructor
  68 }
  69
  70 //_____________________________________________________________________________
  71 AliComparisonDEdx::AliComparisonDEdx(Char_t* name="AliComparisonDEdx", Char_t* title="AliComparisonDEdx",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
  72   AliComparisonObject(name,title),
  73
  74   // dEdx
  75   fDeDxHisto(0),
  76
  77   // Cuts
  78   fCutsRC(0),
  79   fCutsMC(0),
  80
  81   // histogram folder
  82   fAnalysisFolder(0)
  83 {
  84   // named constructor
  85
  86   SetAnalysisMode(analysisMode);
  87   SetHptGenerator(hptGenerator);
  88   Init();
  89 }
  90
  91
  92 //_____________________________________________________________________________
  93 AliComparisonDEdx::~AliComparisonDEdx()
  94 {
  95   // destructor
  96   if(fDeDxHisto)  delete fDeDxHisto; fDeDxHisto=0;
  97
  98   if(fCutsRC) delete fCutsRC; fCutsRC=0;
  99   if(fCutsMC) delete fCutsMC; fCutsMC=0;
 100
 101   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 102 }
 103
 104 //_____________________________________________________________________________
 105 void AliComparisonDEdx::Init()
 106 {
 107   // Init histograms
 108
 109   // TPC dEdx
 110   Int_t nPBins = 31;
 111     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.};
 112     Double_t pMin = 0., pMax = 10.;
 113
 114     if(IsHptGenerator() == kTRUE) {
 115       nPBins = 100;
 116       pMin = 0.; pMax = 100.;
 117     }
 118
 119    //signal:alpha:y:z:snp:tgl:ncls:pid:p
 120    Int_t binsQA[9]    = {600,50, 50,  50, 50, 50, 80, 5, nPBins};
 121    Double_t xminQA[9] = {0,  -4,-20,-250, -1, -2, 0, 0., pMin};
 122    Double_t xmaxQA[9] = {300, 4, 20, 250,  1,  2, 160, 5., pMax};
 123
 124    fDeDxHisto = new THnSparseF("fDeDxHisto","signal:alpha:y:z:snp:tgl:ncls:pid:momentum",9,binsQA,xminQA,xmaxQA);
 125    if(!IsHptGenerator()) fDeDxHisto->SetBinEdges(8,binsP);
 126
 127    fDeDxHisto->GetAxis(0)->SetTitle("signal");
 128    fDeDxHisto->GetAxis(1)->SetTitle("alpha (rad)");
 129    fDeDxHisto->GetAxis(2)->SetTitle("y (cm)");
 130    fDeDxHisto->GetAxis(3)->SetTitle("z (cm)");
 131    fDeDxHisto->GetAxis(4)->SetTitle("snp");
 132    fDeDxHisto->GetAxis(5)->SetTitle("tgl");
 133    fDeDxHisto->GetAxis(6)->SetTitle("ncls");
 134    fDeDxHisto->GetAxis(6)->SetTitle("pid");
 135    fDeDxHisto->GetAxis(7)->SetTitle("p (GeV/c)");
 136    fDeDxHisto->Sumw2();
 137
 138    // Init cuts
 139    if(!fCutsMC)
 140      AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 141    if(!fCutsRC)
 142      AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 143
 144    // init folder
 145    fAnalysisFolder = CreateFolder("folderDEdx","Analysis de/dx Folder");
 146 }
 147
 148 //_____________________________________________________________________________
 149 void AliComparisonDEdx::ProcessTPC(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
 150 {
 151   // Fill dE/dx  comparison information
 152
 153   // Check selection cuts
 154   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 155
 156   Double_t pid = -1;
 157   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetEM() ) pid = 0;
 158   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetMuM() ) pid = 1;
 159   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetPiP() ) pid = 2;
 160   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetKP() ) pid = 3;
 161   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetProt() ) pid = 4;
 162
 163   if (!infoRC->GetESDtrack()) return;
 164   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
 165
 166   Float_t dedx = infoRC->GetESDtrack()->GetTPCsignal();
 167   Int_t ncls = infoRC->GetESDtrack()->GetTPCNcls();
 168
 169   const AliExternalTrackParam *innerParam =  0;
 170   if ((innerParam = infoRC->GetESDtrack()->GetInnerParam()) == 0) return;
 171
 172   Double_t pt = innerParam->Pt();
 173   Double_t lam = TMath::ATan2(innerParam->Pz(),innerParam->Pt());
 174   Double_t p = pt/TMath::Cos(lam);
 175   Double_t alpha = innerParam->GetAlpha();
 176   Double_t y = innerParam->GetY();
 177   Double_t z = innerParam->GetZ();
 178   Double_t snp = innerParam->GetSnp();
 179   Double_t tgl = innerParam->GetTgl();
 180
 181   Double_t vDeDxHisto[9] = {dedx,alpha,y,z,snp,tgl,ncls,pid,p};
 182   fDeDxHisto->Fill(vDeDxHisto);
 183 }
 184
 185 //_____________________________________________________________________________
 186 void AliComparisonDEdx::ProcessTPCITS(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
 187 {
 188   // Fill dE/dx  comparison information
 189
 190    AliDebug(AliLog::kWarning, "Warning: Not implemented");
 191 }
 192
 193 //_____________________________________________________________________________
 194 void AliComparisonDEdx::ProcessConstrained(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
 195 {
 196   // Fill dE/dx  comparison information
 197
 198    AliDebug(AliLog::kWarning, "Warning: Not implemented");
 199 }
 200
 201 //_____________________________________________________________________________
 202 Long64_t AliComparisonDEdx::Merge(TCollection* const list)
 203 {
 204   // Merge list of objects (needed by PROOF)
 205
 206   if (!list)
 207   return 0;
 208
 209   if (list->IsEmpty())
 210   return 1;
 211
 212   TIterator* iter = list->MakeIterator();
 213   TObject* obj = 0;
 214
 215   // collection of generated histograms
 216   Int_t count=0;
 217   while((obj = iter->Next()) != 0)
 218   {
 219     AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);
 220     if (entry == 0) continue;
 221
 222     fDeDxHisto->Add(entry->fDeDxHisto);
 223     count++;
 224   }
 225
 226 return count;
 227 }
 228
 229 //_____________________________________________________________________________
 230 void AliComparisonDEdx::Exec(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
 231 {
 232   // Process comparison information
 233
 234   if(GetAnalysisMode() == 0) ProcessTPC(infoMC,infoRC);
 235   else if(GetAnalysisMode() == 1) ProcessTPCITS(infoMC,infoRC);
 236   else if(GetAnalysisMode() == 2) ProcessConstrained(infoMC,infoRC);
 237   else {
 238     printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 239     return;
 240   }
 241 }
 242
 243 //_____________________________________________________________________________
 244 TH1F* AliComparisonDEdx::MakeResol(TH2F * his, Int_t integ, Bool_t type)
 245 {
 246   // Make resolution histograms
 247   TH1F *hisr, *hism;
 248   if (!gPad) new TCanvas;
 249   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
 250   if (type) return hism;
 251   else
 252     return hisr;
 253 }
 254
 255 //_____________________________________________________________________________
 256 void AliComparisonDEdx::Analyse()
 257 {
 258   // Analyze comparison information and store output histograms
 259   // in the folder "folderDEdx"
 260   //
 261   TH1::AddDirectory(kFALSE);
 262   TObjArray *aFolderObj = new TObjArray;
 263
 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
 309   // export objects to analysis folder
 310   fAnalysisFolder = ExportToFolder(aFolderObj);
 311
 312   // delete only TObjrArray
 313   if(aFolderObj) delete aFolderObj;
 314 }
 315
 316 //_____________________________________________________________________________
 317 TFolder* AliComparisonDEdx::ExportToFolder(TObjArray * array)
 318 {
 319   // recreate folder avery time and export objects to new one
 320   //
 321   AliComparisonDEdx * comp=this;
 322   TFolder *folder = comp->GetAnalysisFolder();
 323
 324   TString name, title;
 325   TFolder *newFolder = 0;
 326   Int_t i = 0;
 327   Int_t size = array->GetSize();
 328
 329   if(folder) {
 330      // get name and title from old folder
 331      name = folder->GetName();
 332      title = folder->GetTitle();
 333
 334          // delete old one
 335      delete folder;
 336
 337          // create new one
 338      newFolder = CreateFolder(name.Data(),title.Data());
 339      newFolder->SetOwner();
 340
 341          // add objects to folder
 342      while(i < size) {
 343            newFolder->Add(array->At(i));
 344            i++;
 345          }
 346   }
 347
 348 return newFolder;
 349 }
 350
 351
 352 //_____________________________________________________________________________
 353 TFolder* AliComparisonDEdx::CreateFolder(TString name,TString title) {
 354 // create folder for analysed histograms
 355 TFolder *folder = 0;
 356   folder = new TFolder(name.Data(),title.Data());
 357
 358   return folder;
 359 }
 360