adding phi, etaphi plots
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 Jun 2008 10:41:36 +0000 (10:41 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 17 Jun 2008 10:41:36 +0000 (10:41 +0000)
adding function to determine acceptance map
adding process type statistics

PWG0/dNdEta/AlidNdEtaCorrectionTask.cxx
PWG0/dNdEta/AlidNdEtaTask.cxx
PWG0/dNdEta/AlidNdEtaTask.h
PWG0/dNdEta/dNdEtaAnalysis.cxx
PWG0/dNdEta/drawPlots.C

index 1dcc6db..9966f14 100644 (file)
@@ -205,10 +205,12 @@ void AlidNdEtaCorrectionTask::CreateOutputObjects()
     fOutput->Add(fDeltaPhi[i]);
   }
 
-  fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 2, -0.5, 1.5, 4, -0.5, 3.5);
+  fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 104, -3.5, 100.5, 4, -0.5, 3.5);
   fOutput->Add(fEventStats);
-  fEventStats->GetXaxis()->SetBinLabel(1, "INEL");
-  fEventStats->GetXaxis()->SetBinLabel(2, "NSD");
+  fEventStats->GetXaxis()->SetBinLabel(1, "INEL"); // x = -3
+  fEventStats->GetXaxis()->SetBinLabel(2, "NSD");  // x = -2
+  for (Int_t i=0; i<100; i++)
+    fEventStats->GetXaxis()->SetBinLabel(4+i, Form("%d", i)); 
 
   fEventStats->GetYaxis()->SetBinLabel(1, "nothing");
   fEventStats->GetYaxis()->SetBinLabel(2, "trg");
@@ -271,7 +273,7 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
 
   if (processType == AliPWG0Helper::kInvalidProcess)
-    AliDebug(AliLog::kError, Form("Unknown process type %d.", processType));
+    AliDebug(AliLog::kError, "Unknown process type.");
 
   // get the MC vertex
   AliGenEventHeader* genHeader = header->GenEventHeader();
@@ -298,10 +300,10 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
   // fill process type
   Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
   // INEL
-  fEventStats->Fill(0.0, biny);
+  fEventStats->Fill(-3, biny);
   // NSD
-  if (processType != AliPWG0Helper::kSD )
-    fEventStats->Fill(1.0, biny);
+  if (processType != AliPWG0Helper::kSD)
+    fEventStats->Fill(-2, biny);
   
   // create list of (label, eta, pt) tuples
   Int_t inputCount = 0;
@@ -446,10 +448,13 @@ void AlidNdEtaCorrectionTask::Exec(Option_t*)
         fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, pt);
 
     // TODO this value might be needed lower for the SPD study (only used for control histograms anyway)
-    if (TMath::Abs(eta) < 1 && pt > 0.2)
+    if (TMath::Abs(eta) < 1.5) // && pt > 0.2)
       nAccepted++;
   }
 
+  if (nAccepted == 0)
+    fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
+
   fMultAll->Fill(nAccepted);
   if (eventTriggered) {
     fMultTr->Fill(nAccepted);
index 73204dd..46e1a02 100644 (file)
@@ -60,6 +60,8 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fdNdEtaAnalysisTracks(0),
   fVertex(0),
   fPartPt(0),
+  fPhi(0),
+  fEtaPhi(0),
   fDeltaPhi(0)
 {
   //
@@ -168,6 +170,15 @@ void AlidNdEtaTask::CreateOutputObjects()
   fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, 10);
   fOutput->Add(fVertexResolution);
 
+  fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
+  fOutput->Add(fPhi);
+
+  fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 80, 0, 2 * TMath::Pi());
+  fOutput->Add(fEtaPhi);
+
+  if (fAnalysisMode == AliPWG0Helper::kSPD)
+    fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 1000, -3.14, 3.14);
+
   if (fReadMC)
   {
     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
@@ -192,9 +203,6 @@ void AlidNdEtaTask::CreateOutputObjects()
     fPartPt->Sumw2();
     fOutput->Add(fPartPt);
   }
-
-  if (fAnalysisMode == AliPWG0Helper::kSPD)
-    fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 1000, -3.14, 3.14);
 }
 
 void AlidNdEtaTask::Exec(Option_t*)
@@ -319,6 +327,12 @@ void AlidNdEtaTask::Exec(Option_t*)
       if (TMath::Abs(deltaPhi) > 1)
         printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
 
+      Float_t phi = mult->GetPhi(i);
+      if (phi < 0)
+        phi += TMath::Pi() * 2;
+      fPhi->Fill(phi);
+      fEtaPhi->Fill(mult->GetEta(i), phi);
+
       fDeltaPhi->Fill(deltaPhi);
 
       etaArr[inputCount] = mult->GetEta(i);
@@ -355,6 +369,12 @@ void AlidNdEtaTask::Exec(Option_t*)
         continue;
       }
 
+      Float_t phi = esdTrack->Phi();
+      if (phi < 0)
+        phi += TMath::Pi() * 2;
+      fPhi->Fill(phi);
+      fEtaPhi->Fill(esdTrack->Eta(), phi);
+
       etaArr[inputCount] = esdTrack->Eta();
       labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
       ptArr[inputCount] = esdTrack->Pt();
@@ -473,7 +493,7 @@ void AlidNdEtaTask::Exec(Option_t*)
       Float_t pt = particle->Pt();
 
        // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison))
-      if (TMath::Abs(eta) < 1.5 && pt > 0.3)
+      if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
         nAcceptedParticles++;
 
       fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
@@ -494,7 +514,7 @@ void AlidNdEtaTask::Exec(Option_t*)
     }
 
     fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
-    if (processType != AliPWG0Helper::kSD )
+    if (processType != AliPWG0Helper::kSD)
       fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
 
     if (eventTriggered)
@@ -552,10 +572,14 @@ void AlidNdEtaTask::Terminate(Option_t *)
   fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
   fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
   fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
-  fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
-  fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
   for (Int_t i=0; i<3; ++i)
     fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
+  fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
+  fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
+
+  fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
+  fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
+  fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
 
   if (!fdNdEtaAnalysisESD)
   {
@@ -629,6 +653,12 @@ void AlidNdEtaTask::Terminate(Option_t *)
   if (fDeltaPhi)
     fDeltaPhi->Write();
 
+  if (fPhi)
+    fPhi->Write();
+
+  if (fEtaPhi)
+    fEtaPhi->Write();
+
   fout->Write();
   fout->Close();
 
@@ -642,6 +672,7 @@ void AlidNdEtaTask::Terminate(Option_t *)
     fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
     fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
     fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
+    fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
 
     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
     {
@@ -666,7 +697,12 @@ void AlidNdEtaTask::Terminate(Option_t *)
     fdNdEtaAnalysisTr->SaveHistograms();
     fdNdEtaAnalysisTrVtx->SaveHistograms();
     fdNdEtaAnalysisTracks->SaveHistograms();
-    fPartPt->Write();
+
+    if (fPartPt)
+      fPartPt->Write();
+
+    if (fVertex)
+      fVertex->Write();
 
     fout->Write();
     fout->Close();
index 8c40dd0..96fd24a 100644 (file)
@@ -10,6 +10,7 @@
 class AliESDtrackCuts;
 class dNdEtaAnalysis;
 class TH1F;
+class TH2F;
 class TH3F;
 class AliESDEvent;
 
@@ -62,9 +63,14 @@ class AlidNdEtaTask : public AliAnalysisTask {
     dNdEtaAnalysis* fdNdEtaAnalysisTr;      //! contains the dndeta from the triggered events
     dNdEtaAnalysis* fdNdEtaAnalysisTrVtx;   //! contains the dndeta from the triggered events with vertex
     dNdEtaAnalysis* fdNdEtaAnalysisTracks;  //! contains the dndeta from the triggered events with vertex counted from the mc particles associated to the tracks (comparing this to the raw values from the esd shows the effect of the detector resolution)
-    // the following are control histograms to check the dNdEtaAnalysis class
+
+    // control histograms (MC)
     TH3F* fVertex;                //! vertex of counted particles
     TH1F* fPartPt;                //! counted particles as function of pt
+
+    // control histograms (ESD)
+    TH1F* fPhi;                   //! raw phi distribution
+    TH2F* fEtaPhi;                //! raw eta - phi distribution
     TH1F* fDeltaPhi;              //! histogram of delta_phi values for tracklets (only for SPD analysis)
 
  private:
index 6614bc2..e0956e3 100644 (file)
@@ -390,29 +390,42 @@ void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, Alid
       Int_t vertexBinBegin = vertexHist->GetXaxis()->FindBin(vertexRangeBegin[vertexPos]);
       Int_t vertexBinEnd   = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]);
 
+      const Int_t *binBegin = 0;
       // adjust acceptance range for SPD
       if (fAnalysisMode == AliPWG0Helper::kSPD)
       {
         //const Int_t binBegin[30] = { 18, 16, 15, 14, 13, 13, 12, 11, 9, 7, 6, 5, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 1 };
-        const Int_t binBegin[30] = { -1, -1, -1, -1, 16, 14, 12, 10, 9, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, -1, -1, -1, -1}; // limit in correction map is 5
+        const Int_t binBeginSPD[30] = { -1, -1, -1, -1, 16, 14, 12, 10, 9, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, -1, -1, -1, -1}; // limit in correction map is 5
         //const Int_t binBegin[30] = { -1, -1, -1, 17, 15, 14, 12, 10, 8, 7, 6, 5, 4, 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, -1, -1, -1};  // limit in correction map is 10
+        binBegin = binBeginSPD;
 
-        Int_t vtxBegin = binBegin[iEta - 1];
-        Int_t vtxEnd   = 18 + 1 - binBegin[30 - iEta];
-
-        // eta range not accessible
-        if (vtxBegin == -1)
-          continue;
-
-        //Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
-        //Printf("%d %d | %d %d", vtxBegin, vertexHist->GetXaxis()->FindBin(GetVtxMin(eta)), vtxEnd, vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta)));
-        //vtxBegin = vertexHist->GetXaxis()->FindBin(GetVtxMin(eta));
-        //vtxEnd = vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta));
-        //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d ", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd);
-        vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin);
-        vertexBinEnd =   TMath::Min(vertexBinEnd, vtxEnd);
-        //Printf("after:  %d %d", vertexBinBegin, vertexBinEnd);
       }
+      else if (fAnalysisMode == AliPWG0Helper::kTPC)
+      {
+//        const Int_t binBegin[30] = { -1, -1, -1, -1, 16, 14, 12, 10, 9, 7, 6, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, -1, -1, -1, -1}; // limit in correction map is 5
+        const Int_t binBeginTPC[30] = {-1, -1, -1, -1, -1, -1, -1, -1, 9, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, -1, -1, -1, -1, -1, -1, -1, -1};
+        binBegin = binBeginTPC;
+      }
+      else if (fAnalysisMode == AliPWG0Helper::kTPCITS)
+      {
+        // TODO create map
+      }
+
+      Int_t vtxBegin = binBegin[iEta - 1];
+      Int_t vtxEnd   = 18 + 1 - binBegin[30 - iEta];
+      
+      // eta range not accessible
+      if (vtxBegin == -1)
+        continue;
+      
+      //Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
+      //Printf("%d %d | %d %d", vtxBegin, vertexHist->GetXaxis()->FindBin(GetVtxMin(eta)), vtxEnd, vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta)));
+      //vtxBegin = vertexHist->GetXaxis()->FindBin(GetVtxMin(eta));
+      //vtxEnd = vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta));
+      //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d ", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd);
+      vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin);
+      vertexBinEnd =   TMath::Min(vertexBinEnd, vtxEnd);
+      //Printf("after:  %d %d", vertexBinBegin, vertexBinEnd);
 
       // no data for this bin
       if (vertexBinBegin > vertexBinEnd)
index 2b944a2..73a28c6 100644 (file)
@@ -1729,7 +1729,7 @@ void DrawTrackletOrigin()
   }
 }
 
-void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
+void DetermineAcceptance(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction", Double_t ptmin=0.2)
 {
   loadlibs();
 
@@ -1739,12 +1739,36 @@ void DetermineAcceptance(const char* fileName = "correction_map.root", const cha
   if (!dNdEtaCorrection->LoadHistograms())
     return;
 
-  hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram();
+  //  hist = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram();
 
-  proj = hist->Project3D("yx");
+  gener = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetGeneratedHistogram();
+  measu = dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram();
+
+  gener->GetZaxis()->SetRange(gener->GetZaxis()->FindBin(ptmin), gener->GetNbinsZ()+1);
+  TH2D *gener_xy = gener->Project3D("yx");
+
+  measu->GetZaxis()->SetRange(measu->GetZaxis()->FindBin(ptmin), measu->GetNbinsZ()+1);
+  TH2D *measu_xy = measu->Project3D("yx");
+
+  cout << measu->GetZaxis()->FindBin(ptmin) << " " << measu->GetNbinsZ()+1 << endl;
+
+  TCanvas *canp = new TCanvas("canp","canp",600,1000);
+  canp->Divide(1,2,0.0001,0.0001);
+  canp->cd(1);
+  gener_xy->Draw("COLZ");
+  canp->cd(2);
+  measu_xy->Draw("COLZ");
+
+
+  TCanvas *canpr = new TCanvas("canpr","canpr",700,500);
+  canpr->cd();
+  TH2D *proj = new TH2D(*gener_xy);
+  proj->Divide(measu_xy);
+
+//   proj = hist->Project3D("yx");
   proj->Draw("COLZ");
 
-  const Float_t limit = 10;
+  const Float_t limit = 5;
 
   TString array = "{";
   TString arrayEnd = "}";