]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG1/AliComparisonDEdx.cxx
QA ref defaut storage setter in sim and rec
[u/mrichter/AliRoot.git] / PWG1 / AliComparisonDEdx.cxx
index effcda226635448a82927593a36cb15657f15f8f..88b2ba42c5c91b95596173641d5828bed497fed7 100644 (file)
@@ -16,7 +16,8 @@
   gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
   LoadMyLibs();
   TFile f("Output.root");
-  AliComparisonDEdx * compObj = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");
+  //AliComparisonDEdx * compObj = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");
+  AliComparisonDEdx * compObj = (AliComparisonDEdx*)cOutput->FindObject("AliComparisonDEdx");
 
   // Analyse comparison data
   compObj->Analyse();
 
 */
 
-#include <iostream>
-
-#include "TFile.h"
-#include "TCint.h"
-#include "TH3F.h"
-#include "TH2F.h"
-#include "TF1.h"
-#include "TProfile.h"
-#include "TProfile2D.h"
-#include "TGraph2D.h"
-#include "TCanvas.h"
-#include "TGraph.h"
-//
+#include <TAxis.h>
+#include <TCanvas.h>
+#include <TH1.h>
+
+#include "AliComparisonDEdx.h" 
 #include "AliESDEvent.h"
-#include "AliESD.h"
-#include "AliESDfriend.h"
-#include "AliESDfriendTrack.h"
-#include "AliRecInfoCuts.h" 
-#include "AliMCInfoCuts.h" 
+#include "AliESDRecInfo.h" 
 #include "AliLog.h" 
-//
+#include "AliMCInfo.h" 
+#include "AliMCInfoCuts.h" 
 #include "AliMathBase.h"
+#include "AliRecInfoCuts.h" 
 #include "AliTreeDraw.h"
-//#include "TStatToolkit.h"
-
-#include "AliMCInfo.h" 
-#include "AliESDRecInfo.h" 
-#include "AliComparisonDEdx.h" 
 
 using namespace std;
 
@@ -71,38 +57,49 @@ AliComparisonDEdx::AliComparisonDEdx():
   AliComparisonObject("AliComparisonDEdx"),
 
   // dEdx 
-  fTPCSignalNormTan(0), 
-  fTPCSignalNormSPhi(0),
-  fTPCSignalNormTPhi(0), 
-  //
-  fTPCSignalNormTanSPhi(0),
-  fTPCSignalNormTanTPhi(0),
-  fTPCSignalNormTanSPt(0), 
+  fDeDxHisto(0),
   
   // Cuts 
   fCutsRC(0), 
   fCutsMC(0),
-  fMCPtMin(0),
-  fMCAbsTanThetaMax(0),
-  fMCPdgCode(0),
 
   // histogram folder 
   fAnalysisFolder(0)
 {
+  // default constructor       
+}
+
+//_____________________________________________________________________________
+AliComparisonDEdx::AliComparisonDEdx(Char_t* name="AliComparisonDEdx", Char_t* title="AliComparisonDEdx",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
+  AliComparisonObject(name,title),
+
+  // dEdx 
+  fDeDxHisto(0),
+  
+  // Cuts 
+  fCutsRC(0), 
+  fCutsMC(0),
+
+  // histogram folder 
+  fAnalysisFolder(0)
+{
+  // named constructor
+
+  SetAnalysisMode(analysisMode);
+  SetHptGenerator(hptGenerator);
   Init();
 }
 
+
 //_____________________________________________________________________________
-AliComparisonDEdx::~AliComparisonDEdx(){
-   
-  if(fTPCSignalNormTan)  delete fTPCSignalNormTan; fTPCSignalNormTan=0; 
-  if(fTPCSignalNormSPhi) delete fTPCSignalNormSPhi; fTPCSignalNormSPhi=0;
-  if(fTPCSignalNormTPhi) delete fTPCSignalNormTPhi; fTPCSignalNormTPhi=0;
-  //
-  if(fTPCSignalNormTanSPhi) delete fTPCSignalNormTanSPhi; fTPCSignalNormTanSPhi=0;
-  if(fTPCSignalNormTanTPhi) delete fTPCSignalNormTanTPhi; fTPCSignalNormTanTPhi=0;
-  if(fTPCSignalNormTanSPt)  delete fTPCSignalNormTanSPt; fTPCSignalNormTanSPt=0;
+AliComparisonDEdx::~AliComparisonDEdx()
+{
+  // destructor
+  if(fDeDxHisto)  delete fDeDxHisto; fDeDxHisto=0; 
 
+  if(fCutsRC) delete fCutsRC; fCutsRC=0;
+  if(fCutsMC) delete fCutsMC; fCutsMC=0;
+  
   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 }
 
@@ -112,92 +109,99 @@ void AliComparisonDEdx::Init()
   // Init histograms
   
   // TPC dEdx
-  fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2,  40,30,70); 
-  fTPCSignalNormTan->SetXTitle("tan(#theta)");
-  fTPCSignalNormTan->SetYTitle("rec. dE/dx / calc. dE/dx");
-
-  fTPCSignalNormSPhi   = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70);
-  fTPCSignalNormSPhi->SetXTitle("sin(#phi)");
-  fTPCSignalNormSPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
-
-  fTPCSignalNormTPhi   = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70); 
-  fTPCSignalNormTPhi->SetXTitle("tan(#phi)");
-  fTPCSignalNormTPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
-
-  fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
-  fTPCSignalNormTanSPhi->SetXTitle("tan(#theta)");
-  fTPCSignalNormTanSPhi->SetYTitle("sin(#phi)");
-  fTPCSignalNormTanSPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
-
-  fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
-  fTPCSignalNormTanTPhi->SetXTitle("tan(#theta)");
-  fTPCSignalNormTanTPhi->SetYTitle("tan(#phi)");
-  fTPCSignalNormTanTPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
-
-  fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3, 40,30,70); 
-  fTPCSignalNormTanSPt->SetXTitle("tan(#theta)");
-  fTPCSignalNormTanSPt->SetYTitle("#sqrt{p_{t}}");
-  fTPCSignalNormTanSPt->SetZTitle("rec. dE/dx / calc. dE/dx");
-
-  // Init cuts
-  if(!fCutsMC) 
-    AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
-  if(!fCutsRC) 
-    AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
-
-    // init folder
-    fAnalysisFolder = CreateFolder("folderDEdx","Analysis de/dx Folder");
+  Int_t nPBins = 31;
+    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.};
+    Double_t pMin = 0., pMax = 10.;
+
+    if(IsHptGenerator() == kTRUE) {
+      nPBins = 100;
+      pMin = 0.; pMax = 100.;
+    }
+
+   //signal:alpha:y:z:snp:tgl:ncls:pid:p
+   Int_t binsQA[9]    = {600,50, 50,  50, 50, 50, 80, 5, nPBins};
+   Double_t xminQA[9] = {0,  -4,-20,-250, -1, -2, 0, 0., pMin};
+   Double_t xmaxQA[9] = {300, 4, 20, 250,  1,  2, 160, 5., pMax};
+
+   fDeDxHisto = new THnSparseF("fDeDxHisto","signal:alpha:y:z:snp:tgl:ncls:pid:momentum",9,binsQA,xminQA,xmaxQA);
+   if(!IsHptGenerator()) fDeDxHisto->SetBinEdges(8,binsP);
+
+   fDeDxHisto->GetAxis(0)->SetTitle("signal");
+   fDeDxHisto->GetAxis(1)->SetTitle("alpha (rad)");
+   fDeDxHisto->GetAxis(2)->SetTitle("y (cm)");
+   fDeDxHisto->GetAxis(3)->SetTitle("z (cm)");
+   fDeDxHisto->GetAxis(4)->SetTitle("snp");
+   fDeDxHisto->GetAxis(5)->SetTitle("tgl");
+   fDeDxHisto->GetAxis(6)->SetTitle("ncls");
+   fDeDxHisto->GetAxis(6)->SetTitle("pid");
+   fDeDxHisto->GetAxis(7)->SetTitle("p (GeV/c)");
+   fDeDxHisto->Sumw2();
+
+   // Init cuts
+   if(!fCutsMC) 
+     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
+   if(!fCutsRC) 
+     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
+
+   // init folder
+   fAnalysisFolder = CreateFolder("folderDEdx","Analysis de/dx Folder");
 }
 
 //_____________________________________________________________________________
-void AliComparisonDEdx::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
-
+void AliComparisonDEdx::ProcessTPC(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
+{
   // Fill dE/dx  comparison information
   
-  Float_t mcpt = infoMC->GetParticle().Pt();
-  Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
-  Float_t mprim = infoMC->GetPrim();
-
-  // distance to Prim. vertex 
-  const Double_t* dv = infoMC->GetVDist(); 
-
-  Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
-  
   // Check selection cuts 
   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
-  if (!isPrim) return;
-  if (infoRC->GetStatus(1)!=3) return;
+  Double_t pid = -1;
+  if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetEM() ) pid = 0;
+  if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetMuM() ) pid = 1;
+  if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetPiP() ) pid = 2;
+  if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetKP() ) pid = 3;
+  if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetProt() ) pid = 4;
+
   if (!infoRC->GetESDtrack()) return;  
   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
-  if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
-  //if (mprim>1.4) return;
-  //if (mprim<0.5) return;
-  if (mprim > fCutsMC->GetMaxTPCSignal()) return;
-  if (mprim < fCutsMC->GetMinTPCSignal()) return;
-  if (infoRC->GetESDtrack()->GetTPCsignalN()<fCutsRC->GetMinTPCsignalN()) return;
-  //
-  Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();
-  Float_t sphi =  infoRC->GetESDtrack()->GetInnerParam()->GetSnp();
-  Float_t tphi =  sphi/TMath::Sqrt((1.-sphi)*(1.+sphi));
-
-  if (TMath::Abs(infoMC->GetParticle().GetPdgCode()) != GetMCPdgCode()) return;
-  //if (mcpt>0.5) {
-  if (mcpt > GetMCPtMin()) {
-    fTPCSignalNormTan->Fill(tantheta,ratio);    // only subset
-  }
 
-  //if (TMath::Abs(tantheta)<0.5){
-  if (TMath::Abs(tantheta) < GetMCAbsTanThetaMax()){
-    fTPCSignalNormSPhi->Fill(sphi,ratio);       // only subset
-    fTPCSignalNormTPhi->Fill(tphi,ratio);       // only subset
-  }
-  fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio);    
-  fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio);    
-  fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);
+  Float_t dedx = infoRC->GetESDtrack()->GetTPCsignal();
+  Int_t ncls = infoRC->GetESDtrack()->GetTPCNcls();
+
+  const AliExternalTrackParam *innerParam =  0;
+  if ((innerParam = infoRC->GetESDtrack()->GetInnerParam()) == 0) return;
+
+  Double_t pt = innerParam->Pt();
+  Double_t lam = TMath::ATan2(innerParam->Pz(),innerParam->Pt());
+  Double_t p = pt/TMath::Cos(lam);
+  Double_t alpha = innerParam->GetAlpha();
+  Double_t y = innerParam->GetY();
+  Double_t z = innerParam->GetZ();
+  Double_t snp = innerParam->GetSnp();
+  Double_t tgl = innerParam->GetTgl();
+
+  Double_t vDeDxHisto[9] = {dedx,alpha,y,z,snp,tgl,ncls,pid,p};
+  fDeDxHisto->Fill(vDeDxHisto); 
 }
 
 //_____________________________________________________________________________
-Long64_t AliComparisonDEdx::Merge(TCollection* list) 
+void AliComparisonDEdx::ProcessTPCITS(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
+{
+  // Fill dE/dx  comparison information
+  
+   AliDebug(AliLog::kWarning, "Warning: Not implemented");
+}
+
+//_____________________________________________________________________________
+void AliComparisonDEdx::ProcessConstrained(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
+{
+  // Fill dE/dx  comparison information
+  
+   AliDebug(AliLog::kWarning, "Warning: Not implemented");
+}
+
+//_____________________________________________________________________________
+Long64_t AliComparisonDEdx::Merge(TCollection* const list) 
 {
   // Merge list of objects (needed by PROOF)
 
@@ -217,14 +221,7 @@ Long64_t AliComparisonDEdx::Merge(TCollection* list)
     AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);
     if (entry == 0) continue;
 
-    fTPCSignalNormTan->Add(entry->fTPCSignalNormTan);
-    fTPCSignalNormSPhi->Add(entry->fTPCSignalNormSPhi);
-    fTPCSignalNormTPhi->Add(entry->fTPCSignalNormTPhi);
-    //
-    fTPCSignalNormTanSPhi->Add(entry->fTPCSignalNormTanSPhi);
-    fTPCSignalNormTanTPhi->Add(entry->fTPCSignalNormTanTPhi);
-    fTPCSignalNormTanSPt->Add(entry->fTPCSignalNormTanSPt);
-
+    fDeDxHisto->Add(entry->fDeDxHisto);
     count++;
   }
 
@@ -232,10 +229,17 @@ return count;
 }
 
 //_____________________________________________________________________________
-void AliComparisonDEdx::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
+void AliComparisonDEdx::Exec(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
 {
   // Process comparison information
-  Process(infoMC,infoRC);
+
+  if(GetAnalysisMode() == 0) ProcessTPC(infoMC,infoRC);
+  else if(GetAnalysisMode() == 1) ProcessTPCITS(infoMC,infoRC);
+  else if(GetAnalysisMode() == 2) ProcessConstrained(infoMC,infoRC);
+  else {
+    printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
+    return;
+  }
 }
 
 //_____________________________________________________________________________
@@ -256,12 +260,10 @@ void AliComparisonDEdx::Analyse()
   // Analyze comparison information and store output histograms
   // in the folder "folderDEdx"
   //
-
   TH1::AddDirectory(kFALSE);
-  
-  AliComparisonDEdx * comp=this;
   TObjArray *aFolderObj = new TObjArray;
 
+  /*
   TH1F *hiss=0;
   TGraph2D * gr=0;
 
@@ -304,6 +306,7 @@ void AliComparisonDEdx::Analyse()
 
   gr->GetHistogram()->SetName("TPCdEdxMeanTanPt_2");
   aFolderObj->Add(gr->GetHistogram());
+  */
 
   // export objects to analysis folder
   fAnalysisFolder = ExportToFolder(aFolderObj);