Add resolution class
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 11 Aug 2008 14:04:37 +0000 (14:04 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 11 Aug 2008 14:04:37 +0000 (14:04 +0000)
TRD/TRDqaRecLinkDef.h
TRD/libTRDqaRec.pkg
TRD/qaRec/AliTRDtrackingResolution.cxx [new file with mode: 0644]
TRD/qaRec/AliTRDtrackingResolution.h [new file with mode: 0644]
TRD/qaRec/run.C

index 2d63f4b..a3282a3 100644 (file)
@@ -7,6 +7,7 @@
 #pragma link off all functions;
 
 #pragma link C++ class  AliTRDtrackInfoGen+;
+#pragma link C++ class  AliTRDtrackingResolution+;
 #pragma link C++ class  AliTRDtrackingEfficiency+;
 #pragma link C++ class  AliTRDtrackingEfficiencyCombined+;
 #pragma link C++ class  AliTRDtrackInfo+;
index b59f988..1a1cf05 100644 (file)
@@ -1,6 +1,7 @@
 #-*- Mode: Makefile -*-
 
 SRCS= qaRec/AliTRDtrackInfoGen.cxx \
+      qaRec/AliTRDtrackingResolution.cxx \
       qaRec/AliTRDtrackingEfficiency.cxx \
       qaRec/AliTRDtrackingEfficiencyCombined.cxx \
       qaRec/AliTRDtrackInfo/AliTRDtrackInfo.cxx
diff --git a/TRD/qaRec/AliTRDtrackingResolution.cxx b/TRD/qaRec/AliTRDtrackingResolution.cxx
new file mode 100644 (file)
index 0000000..f1185d8
--- /dev/null
@@ -0,0 +1,220 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: AliTRDtrackingResolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  Reconstruction QA                                                     //
+//                                                                        //
+//  Authors:                                                              //
+//    Markus Fasel <M.Fasel@gsi.de>                                       //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#include <TObjArray.h>
+#include <TList.h>
+#include <TH1F.h>
+#include <TProfile.h>
+#include <TMath.h>
+
+#include "AliAnalysisManager.h"
+#include "AliTRDseedV1.h"
+#include "AliTrackReference.h"
+#include "TTreeStream.h"
+
+#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
+#include "AliTRDtrackingResolution.h"
+
+#include <cstring>
+
+ClassImp(AliTRDtrackingResolution)
+
+//________________________________________________________
+AliTRDtrackingResolution::AliTRDtrackingResolution(const char * name):
+  AliAnalysisTask(name, ""),
+  fTrackObjects(0x0),
+  fOutputHistograms(0x0),
+  fYres(0x0),
+/*     fZres(0),
+  fYresAngle(0),
+  fPhiRes(0x0),
+  fPhiResAngle(0x0),*/
+  fDebugLevel(0),
+  fDebugStream(0x0)
+{
+//     memset(fYresLayer, 0, sizeof(TH1F *) * kNLayers);
+//     memset(fZresLayer, 0, sizeof(TH1F *) * kNLayers);
+//     memset(fYresLayerAngle, 0, sizeof(TProfile *) * kNLayers);
+//     memset(fPhiResLayer, 0, sizeof(TH1F *) * kNLayers);
+//     memset(fPhiResLayerAngle, 0, sizeof(TProfile *) * kNLayers);
+
+  DefineInput(0, TObjArray::Class());
+  DefineOutput(0, TList::Class());
+}
+
+//________________________________________________________
+void AliTRDtrackingResolution::ConnectInputData(Option_t *){
+  fTrackObjects = dynamic_cast<TObjArray *>(GetInputData(0));
+}
+
+//________________________________________________________
+void AliTRDtrackingResolution::CreateOutputObjects()
+{
+  // spatial resolution
+  printf("Creating Histograms\n");
+  OpenFile(0, "RECREATE");
+  fOutputHistograms = new TList();
+  fYres = new TH1F("fYres", "y-Resolution", 100, -1.5, 1.5);
+  fOutputHistograms->Add(fYres);
+//     fZres = new TH1F("fZres", "z-Resolution", 100, -1.5, 1.5);
+//     fOutputHistograms->Add(fZres);
+// 
+//     fYresAngle = new TProfile("fYresAngle", "y-Resolution - Angluar dependence", 80, -40, 40);
+//     fOutputHistograms->Add(fYresAngle);
+// 
+//     // angular resolution
+//     fPhiRes = new TH1F("fPhiRes", "phi-resolution", 20, -10, 10);
+//     fOutputHistograms->Add(fPhiRes);
+//     
+//     fPhiResAngle = new TProfile("fPhiResAngle", "phi-resolution - Angular dependence", 80, -40, 40);
+//     fOutputHistograms->Add(fPhiResAngle);
+
+/*     for(Int_t iplane = 0; iplane < kNLayers; iplane++){
+    // spatial resolution
+    fYresLayer[iplane] = new TH1F(Form("fYresLayer%d", iplane), Form("y-Resolution in Layer %d", iplane), 100, -1.5, 1.5);
+    fOutputHistograms->Add(fYresLayer[iplane]);
+
+    fZresLayer[iplane] = new TH1F(Form("fZresLayer%d", iplane), Form("z-Resolution in Layer %d", iplane), 100, -1.5, 1.5);
+    fOutputHistograms->Add(fZresLayer[iplane]);
+
+    fYresLayerAngle[iplane] = new TProfile(Form("fYresLayerAngle%d", iplane), Form("y-Resolution in Layer %d - Angluar dependence", iplane), 80, -40, 40);
+    fOutputHistograms->Add(fYresLayerAngle[iplane]);
+
+    // angular resolution
+    fPhiResLayer[iplane] = new TH1F(Form("fPhiResLayer%d", iplane), Form("phi-resolution in Layer %d", iplane), 20, -10, 10);
+    fOutputHistograms->Add(fPhiResLayer[iplane]);
+
+    fPhiResLayerAngle[iplane] = new TProfile(Form("fPhiResAngle%d", iplane), Form("phi-resolution in Layer %d - Angular dependence", iplane), 80, -40, 40);
+    fOutputHistograms->Add(fPhiResLayerAngle[iplane]);
+  }*/
+}
+
+//________________________________________________________
+void AliTRDtrackingResolution::Exec(Option_t *){
+  // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
+  // angular Resolution: res = Tracklet angle - TrackRef Angle
+  Int_t nTrackInfos = fTrackObjects->GetEntriesFast();
+  if(fDebugLevel>=2) printf("Number of Histograms: %d\n", fOutputHistograms->GetEntries());
+  AliTRDtrackInfo *fInfo = 0x0;
+  if(fDebugLevel>=2) printf("Number of TrackInfos: %d\n", nTrackInfos);
+  for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
+    if(fDebugLevel>=2) printf("Doing Object %d\n", iTI);
+    fInfo = dynamic_cast<AliTRDtrackInfo *>(fTrackObjects->UncheckedAt(iTI));
+    // check if ESD and MC-Information are available
+    if(!fInfo || !fInfo->GetTRDtrack() || fInfo->GetNTrackRefs() < 2) continue; 
+    AliTRDseedV1*fTracklet = 0;
+    AliTrackReference *fTrackRefs[2];
+    for(Int_t iplane = 0; iplane < kNLayers; iplane++){
+      if(fDebugLevel>=2) printf("plane %d\n", iplane);
+      fTracklet = fInfo->GetTracklet(iplane);
+      if(!fTracklet) continue;
+      // check for 2 track ref where the radial position has a distance less than 3.7mm
+      if(fDebugLevel>=2) printf("Find Track References for x = %f\n", fTracklet->GetX0());
+      if(fDebugLevel>=2) printf("Number of Clusters: %d\n", fTracklet->GetN());
+      Int_t nFound = 0;
+      memset(fTrackRefs, 0, sizeof(AliTrackReference*) * 2);
+      AliTrackReference *tempTrackRef = 0;
+      for(Int_t itr = 0; itr < fInfo->GetNTrackRefs(); itr++){
+        if(fDebugLevel>=2) printf("nFound = %d\n", nFound);
+        if(nFound >= 2) break;
+        tempTrackRef = fInfo->GetTrackRef(itr);
+        if(!tempTrackRef) continue;
+        if(fDebugLevel>=2) printf("TrackRef %d: x = %f\n", itr, tempTrackRef->LocalX());
+        if(fTracklet->GetX0() - tempTrackRef->LocalX() > 3.7) continue;
+        if(tempTrackRef->LocalX() - fTracklet->GetX0() > 3.7) break;
+        if(fDebugLevel>=2) printf("accepted\n");
+        if(nFound == 1)
+          if(fTrackRefs[0]->LocalX() >= tempTrackRef->LocalX()) continue;
+        fTrackRefs[nFound++] = tempTrackRef;
+      }
+      if(fDebugLevel>=2) printf("nFound = %d\n", nFound);
+      if(nFound < 2) continue;
+      // We found 2 track refs for the tracklet, get y and z at the anode wire by a linear approximation
+      Double_t dx = fTrackRefs[1]->LocalX() - fTrackRefs[0]->LocalX();
+      Double_t dydx = (fTrackRefs[1]->LocalY() - fTrackRefs[0]->LocalY()) / dx;
+      Double_t dzdx = (fTrackRefs[1]->Z() - fTrackRefs[0]->Z()) / dx;
+      Double_t dx0 = fTrackRefs[1]->LocalX() - fTracklet->GetX0();
+      Double_t ymc =  fTrackRefs[1]->LocalY() - dydx*dx0;
+      Double_t zmc =  fTrackRefs[1]->Z() - dzdx*dx0;
+      
+      Double_t dy = fTracklet->GetYfit(0) - ymc;
+      Double_t dz = fTracklet->GetZfit(0) - zmc;
+      //res_y *= 100; // in mm
+      Double_t momentum = fTrackRefs[0]->P();
+      
+      // Fill Histograms
+      if(fDebugLevel>=2) printf("dy = %f\n", dy);
+      fYres->Fill(dy);
+//                     fZres->Fill(res_z);
+/*                     fYresLayer[iplane]->Fill(res_y);
+      fZresLayer[iplane]->Fill(res_z);*/
+
+//       Double_t phi     = fTrackRefs[0]->Phi();
+//       Double_t theta   = fTrackRefs[0]->Theta();
+      Double_t phi   = TMath::ATan(dydx);
+      Double_t theta = TMath::ATan(dzdx);
+      
+//                     fYresAngle->Fill(phi, res_y);
+//                     fYresLayerAngle[iplane]->Fill(phi, res_y);
+      
+      Double_t dphi   = TMath::ATan(fTracklet->GetZfit(1)) - phi;
+      
+//                     fPhiRes->Fill(dphi);
+//                     fPhiResLayer[iplane]->Fill(dphi);
+//                     fPhiResAngle->Fill(phi, dphi);
+//                     fPhiResLayerAngle[iplane]->Fill(phi, dphi);
+      
+      // Fill Debug Tree
+      if(fDebugLevel>=1){
+        (*fDebugStream) << "Resolution"
+          << "plane="          << iplane
+          << "p="       << momentum
+          << "dx="      << dx
+          << "dy="               << dy
+          << "dz="               << dz
+          << "phi="                    << phi
+          << "theta="          << theta
+          << "dphi="           << dphi
+          << "\n";
+      }
+    }
+  }
+  PostData(0, fOutputHistograms);
+}
+
+//________________________________________________________
+void AliTRDtrackingResolution::Terminate(Option_t *){
+  //printf("Tracking Resolution: Terminate\n");
+  if(fDebugStream) delete fDebugStream;
+}
+
+//________________________________________________________
+void AliTRDtrackingResolution::SetDebugLevel(Int_t level){
+  fDebugLevel = level;
+  if(!fDebugLevel) return;
+  if(fDebugStream) return;
+  fDebugStream = new TTreeSRedirector("TRD.Resolution.root");
+}
diff --git a/TRD/qaRec/AliTRDtrackingResolution.h b/TRD/qaRec/AliTRDtrackingResolution.h
new file mode 100644 (file)
index 0000000..eb77816
--- /dev/null
@@ -0,0 +1,60 @@
+#ifndef ALITRDTRACKINGRESOLUTION_H
+#define ALITRDTRACKINGRESOLUTION_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliTRDtrackingResolution.h 27496 2008-07-22 08:35:45Z cblume $ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  Reconstruction QA                                                     //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#include "AliAnalysisTask.h"
+
+class TObjArray;
+class TList;
+class TH1F;
+class TProfile;
+class TTreeSRedirector;
+
+class AliTRDtrackingResolution : public AliAnalysisTask{
+public:
+  AliTRDtrackingResolution(const char *name = "TRD Tracking Resolution");
+  ~AliTRDtrackingResolution(){};
+  
+  void ConnectInputData(Option_t *);
+  void CreateOutputObjects();
+  void Exec(Option_t *);
+  void Terminate(Option_t *);
+  
+  void SetDebugLevel(Int_t level);
+  Int_t GetDebugLevel() const { return fDebugLevel;}
+  
+private:
+  AliTRDtrackingResolution(const AliTRDtrackingResolution&);
+  AliTRDtrackingResolution& operator=(const AliTRDtrackingResolution&);
+
+private:
+  enum{
+    kNLayers = 6
+  };
+  TObjArray *fTrackObjects;     // Input Track Info container
+  TList     *fOutputHistograms; // Container for the output histograms
+  TH1F      *fYres;             // y-Resolution
+/*             TH1F *fZres;                                                                                                            // z-Resolution
+    TProfile *fYresAngle;                                                                              // y-Resolution dependent on Angle*/
+/*             TH1F *fYresLayer[kNLayers];                                                     // y-Resolution histograms for each Layer
+    TH1F *fZresLayer[kNLayers];                                                        // z-Resolution histograms for each Layer
+    TProfile *fYresLayerAngle[kNLayers];               // y-Resolution histograms for each Layer - Angular Dependence*/
+/*             TH1F *fPhiRes;                                                                                                  // Angular resolution in Phi-Direction
+    TProfile *fPhiResAngle;                                                                    // Phi-resolution dependent on angle*/
+/*             TH1F *fPhiResLayer[kNLayers];                                           // Phi-Resolution histograms for each Layer
+    TProfile *fPhiResLayerAngle[kNLayers];     // Phi-resolution histograms for each Layer - Angular Dependence */
+    Int_t fDebugLevel;                                                                                 // Debug Level
+    TTreeSRedirector *fDebugStream;                            // Debug stream
+  
+  ClassDef(AliTRDtrackingResolution, 1) // tracking resolution task
+};
+#endif
index a385231..7f62c43 100644 (file)
@@ -1,15 +1,68 @@
-void run(const Char_t *list=0x0, Int_t nmax=-1) {
+// Steer TRD QA train for Reconstruction (Clusterizer, Tracking and PID).
+// 
+// Usage:
+//   run.C(tasks, files, entries)
+//   tasks : "ALL" or one/more of the following:
+//     "EFF"  : TRD Tracking Efficiency 
+//     "EFFC" : TRD Tracking Efficiency Combined (barrel + stand alone) - only in case of simulations
+//     "RES"  : TRD tracking Resolution
+//     "PID"  : TRD PID - pion efficiency 
+//     "PIDR" : TRD PID - reference data
+// 
+// Authors:
+//   Alex Bercuci (A.Bercuci@gsi.de) 
+//   Markus Fasel (m.Fasel@gsi.de) 
+
+#define BIT(n)       (1 << (n))
+#define SETBIT(n,i)  ((n) |= BIT(i))
+#define TESTBIT(n,i) ((Bool_t)(((n) & BIT(i)) != 0))
+
+const Int_t fknTasks = 3;
+Char_t *fTaskName[fknTasks] = {"Barrel Tracking Effiency", "Combined Tracking Efficiency", "Tracking Resolution"};
+enum AliTRDrecoTasks{
+  kTrackingEfficiency = 0
+  ,kTrackingCombinedEfficiency = 1
+  ,kTrackingResolution = 2
+};
+void run(Char_t *tasks="ALL", const Char_t *files=0x0, Int_t nmax=-1)
+{
+
+
   TStopwatch timer;
   timer.Start();
 
   gSystem->Load("libANALYSIS.so");
   gSystem->Load("libTRDqaRec.so");
-       
+  
+  Int_t fSteerTask = 0; 
+  TObjArray *task = TString(tasks).Tokenize(" ");
+  for(Int_t isel = 0; isel < task->GetEntriesFast(); isel++){
+    TString s = (dynamic_cast<TObjString *>(task->UncheckedAt(isel)))->String();
+    if(s.CompareTo("ALL") == 0){
+      for(Int_t itask = 0; itask < fknTasks; itask++) SETBIT(fSteerTask, itask);
+      continue;
+    } else if(s.CompareTo("EFF") == 0){
+      SETBIT(fSteerTask, kTrackingEfficiency);
+      continue;
+    } else if(s.CompareTo("EFFC") == 0){
+      SETBIT(fSteerTask, kTrackingCombinedEfficiency);
+      continue;
+    } else if(s.CompareTo("RES" ) == 0){
+      SETBIT(fSteerTask, kTrackingResolution);
+      continue;
+    } else{
+      Info("run.C", Form("Task %s not implemented (yet).", s.Data()));
+      continue;
+    }
+  }
+  printf("\n\tRUNNING TRAIN FOR TASKS:\n");
+  for(itask = 0; itask < fknTasks; itask++){
+    if(TESTBIT(fSteerTask, itask)) printf("\t%s\n", fTaskName[itask]);
+  }
 
   //____________________________________________//
-
   gROOT->LoadMacro(Form("%s/TRD/qaRec/CreateESDChain.C", gSystem->ExpandPathName("$ALICE_ROOT")));
-  TChain *chain = CreateESDChain(list, nmax);
+  TChain *chain = CreateESDChain(files, nmax);
   //chain->SetBranchStatus("*", 0);
   chain->SetBranchStatus("*FMD*",0);
   chain->SetBranchStatus("*Calo*",0);
@@ -44,25 +97,39 @@ void run(const Char_t *list=0x0, Int_t nmax=-1) {
 
   //____________________________________________
   // TRD barrel tracking efficiency
-  AliTRDtrackingEfficiency *task2 = new AliTRDtrackingEfficiency();
-  task2->SetDebugLevel(1);
-  mgr->AddTask(task2);
-  //Create containers for input/output
-  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("TrackingEfficiency", TList::Class(), AliAnalysisManager::kOutputContainer, "TRD.TrackingEfficiency.root");
-  mgr->ConnectInput( task2, 0, coutput1);
-  mgr->ConnectOutput(task2, 0, coutput2);
+  if(TESTBIT(fSteerTask, kTrackingEfficiency)){
+    AliTRDtrackingEfficiency *task2 = new AliTRDtrackingEfficiency();
+    task2->SetDebugLevel(1);
+    mgr->AddTask(task2);
+    //Create containers for input/output
+    AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("TrackingEfficiency", TList::Class(), AliAnalysisManager::kOutputContainer, "TRD.TrackingEfficiency.root");
+    mgr->ConnectInput( task2, 0, coutput1);
+    mgr->ConnectOutput(task2, 0, coutput2);
+  }
 
   //____________________________________________
   // TRD combined tracking efficiency
-  AliTRDtrackingEfficiencyCombined *task3 = new AliTRDtrackingEfficiencyCombined();
-  task3->SetDebugLevel(0);
-  mgr->AddTask(task3);
-  // Create containers for input/output
-  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("TrackingEfficiencyCombined", TObjArray::Class(), AliAnalysisManager::kOutputContainer, "TRD.TrackingEfficiencyCombined.root");
-  mgr->ConnectInput( task3, 0, coutput1);
-  mgr->ConnectOutput(task3, 0, coutput3);
-
+  if(TESTBIT(fSteerTask, kTrackingCombinedEfficiency)){
+    AliTRDtrackingEfficiencyCombined *task3 = new AliTRDtrackingEfficiencyCombined();
+    task3->SetDebugLevel(0);
+    mgr->AddTask(task3);
+    // Create containers for input/output
+    AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("TrackingEfficiencyCombined", TObjArray::Class(), AliAnalysisManager::kOutputContainer, "TRD.TrackingEfficiencyCombined.root");
+    mgr->ConnectInput( task3, 0, coutput1);
+    mgr->ConnectOutput(task3, 0, coutput3);
+  }
 
+  //____________________________________________
+  // TRD combined tracking efficiency
+  if(TESTBIT(fSteerTask, kTrackingResolution)){
+    AliTRDtrackingResolution *task4 = new AliTRDtrackingResolution();
+    task4->SetDebugLevel(1);
+    mgr->AddTask(task4);
+    // Create containers for input/output
+    AliAnalysisDataContainer *coutput4 = mgr->CreateContainer("Tracking Resolution", TList::Class(), AliAnalysisManager::kOutputContainer, "TRD.TrackingResolution.root");
+    mgr->ConnectInput( task4, 0, coutput1);
+    mgr->ConnectOutput(task4, 0, coutput4);
+  }
 
   if (!mgr->InitAnalysis()) return;
   mgr->PrintStatus();