updates in the cluster error parameterization task

diff --git a/TRD/qaRec/AliTRDclusterResolution.cxx b/TRD/qaRec/AliTRDclusterResolution.cxx
index 85e71573009268a49b73df2df303030d1f26dcee..5d93b51dc711a3f12aac74a6e244e1988db80595 100644 (file)
--- a/TRD/qaRec/AliTRDclusterResolution.cxx
+++ b/TRD/qaRec/AliTRDclusterResolution.cxx
@@ -1,5 +1,7 @@
 #include "AliTRDclusterResolution.h"
 #include "AliTRDtrackInfo/AliTRDclusterInfo.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDpadPlane.h"
 
 #include "AliLog.h"
 
@@ -15,7 +17,7 @@ ClassImp(AliTRDclusterResolution)
 
 //_______________________________________________________
 AliTRDclusterResolution::AliTRDclusterResolution()
-  : AliTRDrecoTask("CalibClRes", "Calibrate Cluster Resolution for Tracking")
+  : AliTRDrecoTask("ClErrParam", "Cluster Error Parametrization")
   ,fInfo(0x0)
   ,fResults(0x0)
   ,fAt(0x0)
@@ -94,17 +96,38 @@ void AliTRDclusterResolution::GetRefFigure(Int_t ifig)
 TObjArray* AliTRDclusterResolution::Histos()
 {
   if(fContainer) return fContainer;
-  fContainer = new TObjArray(kN+1);
+  fContainer = new TObjArray(3);
   //fContainer->SetOwner(kTRUE);
 
+  TH2I *h2 = 0x0;
+  TObjArray *arr = 0x0;
+
+  fContainer->AddAt(h2 = new TH2I("h_q", "", 50, 2.2, 7.5, 100, -.5, .5), kQRes);
+  h2->SetXTitle("log(q) [a.u.]");
+  h2->SetYTitle("#Delta y[cm]");
+  h2->SetZTitle("entries");
+
+  Double_t w = 0.;
+  AliTRDgeometry geo;
+  fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kYRes);
+  for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
+    w = .5*geo.GetPadPlane(ily, 2)->GetWidthIPad();
+    arr->AddAt(h2 = new TH2I(Form("h_y%d", ily), Form("Ly[%d]", ily), 50, -w, w, 100, -.5, .5), ily);
+    h2->SetXTitle("y_{w} [w]");
+    h2->SetYTitle("#Delta y[cm]");
+    h2->SetZTitle("entries");
+  }
+
+  fContainer->AddAt(arr = new TObjArray(kN), kSXRes);
   Int_t ih = 0;
   for(Int_t id=1; id<=fAd->GetNbins(); id++){
     for(Int_t it=1; it<=fAt->GetNbins(); it++){
-      fContainer->AddAt(new TH2I(Form("h_d%02dt%02d", id, it), Form("d_{wire}(%5.2f-%5.2f)[mm] x_{drift}(%5.2f-%5.2f)[mm]", 10.*fAd->GetBinCenter(id)- 5.*fAd->GetBinWidth(id), 10.*fAd->GetBinCenter(id)+ 5.*fAd->GetBinWidth(id), 10.*fAt->GetBinCenter(it)- 5.*fAt->GetBinWidth(it), 10.*fAt->GetBinCenter(it)+ 5.*fAt->GetBinWidth(it)), 30, -.15, .15, 100, -.5, .5), ih++);
+      arr->AddAt(h2 = new TH2I(Form("h_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] x_{drift}(%5.2f-%5.2f)[mm]", 10.*fAd->GetBinCenter(id)- 5.*fAd->GetBinWidth(id), 10.*fAd->GetBinCenter(id)+ 5.*fAd->GetBinWidth(id), 10.*fAt->GetBinCenter(it)- 5.*fAt->GetBinWidth(it), 10.*fAt->GetBinCenter(it)+ 5.*fAt->GetBinWidth(it)), 30, -.15, .15, 100, -.5, .5), ih++);
+      h2->SetXTitle("tg#phi");
+      h2->SetYTitle("#Delta y[cm]");
+      h2->SetZTitle("entries");
     }
   }
-  fContainer->AddAt(new TH2I("h_q", "", 50, 2.2, 7.5, 100, -.5, .5), ih++);
-  fContainer->AddAt(new TH2I("h_y", "", 50, -.35, .35, 100, -.5, .5), ih++);
   return fContainer;
 }
 
@@ -114,6 +137,10 @@ void AliTRDclusterResolution::Exec(Option_t *)
   Int_t det, t;
   Float_t x, y, z, q, dy, dydx, dzdx, cov[3], covcl[3];
   TH2I *h2 = 0x0;
+
+  TObjArray *arr0 = (TObjArray*)fContainer->At(kYRes);
+  TObjArray *arr1 = (TObjArray*)fContainer->At(kSXRes);
+
   const AliTRDclusterInfo *cli = 0x0;
   TIterator *iter=fInfo->MakeIterator();
   while((cli=dynamic_cast<AliTRDclusterInfo*>((*iter)()))){
@@ -128,7 +155,7 @@ void AliTRDclusterResolution::Exec(Option_t *)
       AliWarning(Form("Distance to anode %f outside allowed range", cli->GetAnisochronity()));
       continue;
     }
-    if(!(h2 = (TH2I*)fContainer->At((id-1)*kNTB+it-1))){
+    if(!(h2 = (TH2I*)arr1->At((id-1)*kNTB+it-1))){
       AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", (id-1)*kNTB+it-1, id, it, cli->GetDriftLength(), cli->GetAnisochronity()));
       continue;
     }
@@ -142,10 +169,10 @@ void AliTRDclusterResolution::Exec(Option_t *)
     // only for phi equal exB 
     if(TMath::Abs(dydx-fExB)<.01){
       cli->GetCluster(det, x, y, z, q, t, covcl);
-      h2 = (TH2I*)fContainer->At(kN);
+      h2 = (TH2I*)fContainer->At(kQRes);
       h2->Fill(TMath::Log(q), dy);
       
-      h2 = (TH2I*)fContainer->At(kN+1);
+      h2 = (TH2I*)arr0->At(AliTRDgeometry::GetLayer(det));
       h2->Fill(cli->GetYDisplacement(), dy);
     }
   }
@@ -158,12 +185,12 @@ Bool_t AliTRDclusterResolution::PostProcess()
 {
   if(!fContainer) return kFALSE;
   
+  TObjArray *arr = 0x0;
   TH2 *h2 = 0x0; 
   if(!fResults){
     TGraphErrors *g = 0x0;
     fResults = new TObjArray(4);
     fResults->SetOwner();
-    TObjArray *arr = 0x0;
     fResults->AddAt(arr = new TObjArray(2), kQRes);
     arr->SetOwner();
     arr->AddAt(g = new TGraphErrors(), 0);
@@ -193,7 +220,7 @@ Bool_t AliTRDclusterResolution::PostProcess()
   } else {
     TObject *o = 0x0;
     TIterator *iter=fResults->MakeIterator();
-    while((o=((*iter)()))) o->Clear();
+    while((o=((*iter)()))) o->Clear(); // maybe it is wrong but we should never reach this point
   }
 
   TF1 f("f", "gaus", -.5, .5);
@@ -203,7 +230,7 @@ Bool_t AliTRDclusterResolution::PostProcess()
   TH1D *h1 = 0x0;
 
   // process resolution dependency on charge
-  if((h2 = (TH2I*)fContainer->At(kN))) {
+  if((h2 = (TH2I*)fContainer->At(kQRes))) {
     TObjArray *arr = (TObjArray*)fResults->At(kQRes);
     TGraphErrors *gqm = (TGraphErrors*)arr->At(0);
     TGraphErrors *gqs = (TGraphErrors*)arr->At(1);
@@ -225,71 +252,76 @@ Bool_t AliTRDclusterResolution::PostProcess()
   } else AliWarning("Missing dy=f(Q) histo");
 
   // process resolution dependency on y displacement
-  if((h2 = (TH2I*)fContainer->At(kN+1))) {
-    TObjArray *arr = (TObjArray*)fResults->At(kYRes);
-    TGraphErrors *gym = (TGraphErrors*)arr->At(0);
-    TGraphErrors *gys = (TGraphErrors*)arr->At(1);
-    ax = h2->GetXaxis();
-    for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
-      Float_t yd = ax->GetBinCenter(ix);
-      h1 = h2->ProjectionY("py", ix, ix);
-      if(h1->GetEntries() < 50) continue;
-      Adjust(&f, h1);
-      h1->Fit(&f, "Q");
-  
-      // Fill sy^2 = f(log(q))
-      Int_t ip = gym->GetN();
-      gym->SetPoint(ip, yd, 10.*f.GetParameter(1));
-      gym->SetPointError(ip, 0., 10.*f.GetParError(1));
-      gys->SetPoint(ip, yd, 1.e2*f.GetParameter(2)*f.GetParameter(2));
-      gys->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2));
-    } 
-  } else AliWarning("Missing dy=f(y_d) histo");
-
-  // process resolution dependency on drift legth and drift cell width
-  TGraphErrors *gm = new TGraphErrors(), *gs = new TGraphErrors();
-  Float_t d(0.), x(0.);
-  TH2F *hsx = (TH2F*)fResults->At(kSXRes);
-  TH2F *hsy = (TH2F*)fResults->At(kSYRes);
-  for(Int_t id=1; id<=fAd->GetNbins(); id++){
-    d = fAd->GetBinCenter(id); //[mm]
-    printf("Doing d=%f[mm]\n", d);
-    for(Int_t it=1; it<=fAt->GetNbins(); it++){
-      x = fAt->GetBinCenter(it); //[mm]
-      printf("Doing x=%f[mm]\n", x);
-      Int_t idx = (id-1)*kNTB+it-1;
-      if(!(h2 = (TH2I*)fContainer->At(idx))) {
-        AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", idx, id, it, x, d));
-        continue;
-      }
-
-      Int_t ip = 0;
+  if((arr = (TObjArray*)fContainer->At(kYRes))) {
+    for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
+      if(!(h2 = (TH2I*)arr->At(ily))) continue;
+      TObjArray *arrg = (TObjArray*)fResults->At(kYRes);
+      TGraphErrors *gym = (TGraphErrors*)arrg->At(0);
+      TGraphErrors *gys = (TGraphErrors*)arrg->At(1);
       ax = h2->GetXaxis();
       for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
-        Float_t dydx = ax->GetBinCenter(ix) - fExB;
-        if(dydx<0.) continue;
+        Float_t yd = ax->GetBinCenter(ix);
         h1 = h2->ProjectionY("py", ix, ix);
         if(h1->GetEntries() < 50) continue;
         Adjust(&f, h1);
-        h1->Fit(&f, "Q", "", -.2, .2);
+        h1->Fit(&f, "Q");
+    
+        // Fill sy^2 = f(log(q))
+        Int_t ip = gym->GetN();
+        gym->SetPoint(ip, yd, 10.*f.GetParameter(1));
+        gym->SetPointError(ip, 0., 10.*f.GetParError(1));
+        gys->SetPoint(ip, yd, 1.e2*f.GetParameter(2)*f.GetParameter(2));
+        gys->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2));
+      } 
+    }
+  } else AliWarning("Missing dy=f(y_d) container");
 
-        // Fill sy^2 = f(tg_phi^2)
-        gm->SetPoint(ip, dydx, 10.*f.GetParameter(1));
-        gm->SetPointError(ip, 0., 10.*f.GetParError(1));
-        gs->SetPoint(ip, dydx*dydx, 1.e2*f.GetParameter(2)*f.GetParameter(2));
-        gs->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2));
-        ip++;
-      }
-      for(;ip<gm->GetN(); ip++){
-        gm->RemovePoint(ip);gs->RemovePoint(ip);
+  // process resolution dependency on drift legth and drift cell width
+  TGraphErrors *gm = new TGraphErrors(), *gs = new TGraphErrors();
+  Float_t d(0.), x(0.);
+  TH2F *hsx = (TH2F*)fResults->At(kSXRes);
+  TH2F *hsy = (TH2F*)fResults->At(kSYRes);
+  if((arr = (TObjArray*)fContainer->At(kSXRes))) {
+    for(Int_t id=1; id<=fAd->GetNbins(); id++){
+      d = fAd->GetBinCenter(id); //[mm]
+      printf("Doing d=%f[mm]\n", d);
+      for(Int_t it=1; it<=fAt->GetNbins(); it++){
+        x = fAt->GetBinCenter(it); //[mm]
+        printf("Doing x=%f[mm]\n", x);
+        Int_t idx = (id-1)*kNTB+it-1;
+        if(!(h2 = (TH2I*)arr->At(idx))) {
+          AliWarning(Form("Missing histo at index idx[%3d] [id[%2d] it[%2d]] xd[%f] d[%f]\n", idx, id, it, x, d));
+          continue;
+        }
+  
+        Int_t ip = 0;
+        ax = h2->GetXaxis();
+        for(Int_t ix=1; ix<=ax->GetNbins(); ix++){
+          Float_t dydx = ax->GetBinCenter(ix) - fExB;
+          if(dydx<0.) continue;
+          h1 = h2->ProjectionY("py", ix, ix);
+          if(h1->GetEntries() < 50) continue;
+          Adjust(&f, h1);
+          h1->Fit(&f, "Q", "", -.2, .2);
+  
+          // Fill sy^2 = f(tg_phi^2)
+          gm->SetPoint(ip, dydx, 10.*f.GetParameter(1));
+          gm->SetPointError(ip, 0., 10.*f.GetParError(1));
+          gs->SetPoint(ip, dydx*dydx, 1.e2*f.GetParameter(2)*f.GetParameter(2));
+          gs->SetPointError(ip, 0., 2.e2*f.GetParameter(2)*f.GetParError(2));
+          ip++;
+        }
+        for(;ip<gm->GetN(); ip++){
+          gm->RemovePoint(ip);gs->RemovePoint(ip);
+        }
+        gs->Fit(&sig, "Q");
+        hsx->SetBinContent(id, it, sig.GetParameter(1));
+        hsx->SetBinError(id, it, sig.GetParError(1));
+        hsy->SetBinContent(id, it, sig.GetParameter(0));
+        hsy->SetBinError(id, it, sig.GetParError(0));
       }
-      gs->Fit(&sig, "Q");
-      hsx->SetBinContent(id, it, sig.GetParameter(1));
-      hsx->SetBinError(id, it, sig.GetParError(1));
-      hsy->SetBinContent(id, it, sig.GetParameter(0));
-      hsy->SetBinError(id, it, sig.GetParError(0));
     }
-  }
+  } else AliWarning("Missing dy=f(x_d, d_wire) container");
   delete gm; delete gs;
 
   return kTRUE;

diff --git a/TRD/qaRec/run.C b/TRD/qaRec/run.C
index 50ce720bca527448814b47cc7cd326732af3c7e3..9a9fff1cdb78f5276caf7ab607395e9b28058d6f 100644 (file)
--- a/TRD/qaRec/run.C
+++ b/TRD/qaRec/run.C
@@ -6,6 +6,7 @@
 //     "EFF"  : TRD Tracking Efficiency 
 //     "EFFC" : TRD Tracking Efficiency Combined (barrel + stand alone) - only in case of simulations
 //     "RES"  : TRD tracking Resolution
+//     "CLRES": clusters Resolution
 //     "CAL"  : TRD calibration
 //     "PID"  : TRD PID - pion efficiency 
 //     "PIDR" : TRD PID - reference data
@@ -57,6 +58,7 @@
 #include "TRD/qaRec/AliTRDpidChecker.h"
 #include "TRD/qaRec/AliTRDpidRefMaker.h"
 #include "TRD/qaRec/AliTRDcheckDetector.h"
+#include "TRD/qaRec/AliTRDclusterResolution.h"
 #endif
 
 #include "run.h"
@@ -96,7 +98,7 @@ void run(Char_t *tasks="ALL", const Char_t *files=0x0, Int_t nmax=-1)
       fHasMCdata = kFALSE;
     } else { 
       Bool_t foundOpt = kFALSE;  
-      for(Int_t itask = 1; itask < NTRDTASKS; itask++){
+      for(Int_t itask = 1; itask <= NTRDTASKS; itask++){
         if(s.CompareTo(fgkTRDtaskOpt[itask]) != 0) continue;
         SETBIT(fSteerTask, itask);
         foundOpt = kTRUE;
@@ -105,6 +107,9 @@ void run(Char_t *tasks="ALL", const Char_t *files=0x0, Int_t nmax=-1)
       if(!foundOpt) Info("run.C", Form("Task %s not implemented (yet).", s.Data()));
     }
   }
+  // extra rules for calibration tasks
+  if(TSTBIT(fSteerTask, kClusterErrorParam)) SETBIT(fSteerTask, kTrackingResolution);
+
   // define task list pointers;
   AliTRDrecoTask *taskPtr[NTRDTASKS], *task = 0x0;
   memset(taskPtr, 0, NTRDTASKS*sizeof(AliAnalysisTask*));
@@ -197,19 +202,26 @@ void run(Char_t *tasks="ALL", const Char_t *files=0x0, Int_t nmax=-1)
     mgr->ConnectInput( task, 0, coutput1);
     mgr->ConnectOutput(task, 0, mgr->CreateContainer(task->GetName(), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Task%s.root", task->GetName())));
 
-    AliAnalysisDataContainer *co1 = mgr->CreateContainer(Form("%sClRez", task->GetName()), TObjArray::Class(), AliAnalysisManager::kExchangeContainer);
-    mgr->ConnectOutput(task, 1, co1);
-
-    AliAnalysisDataContainer *co2 = mgr->CreateContainer(Form("%sClRes", task->GetName()), TObjArray::Class(), AliAnalysisManager::kExchangeContainer);
-    mgr->ConnectOutput(task, 2, co2);
-
-    AliAnalysisDataContainer *co3 = mgr->CreateContainer(Form("%sTrkltRes", task->GetName()), TObjArray::Class(), AliAnalysisManager::kExchangeContainer);
-    mgr->ConnectOutput(task, 3, co3);
+    // Create output containers for calibration tasks
+    const Int_t nc = 3;
+    const Char_t *cn[nc] = {"ClRez", "ClRes", "TrkltRes"}; 
+    AliAnalysisDataContainer *co[nc]; 
+    for(Int_t ic = 0; ic<nc; ic++){
+      co[ic] = mgr->CreateContainer(Form("%s%s", task->GetName(), cn[ic]), TObjArray::Class(), AliAnalysisManager::kExchangeContainer);
+      mgr->ConnectOutput(task, 1+ic, co[ic]);
+    }
     
     // test reconstruction calibration plugin
-    mgr->AddTask(task = new AliTRDclusterResolution());
-    mgr->ConnectInput(task, 0, co2);
-    mgr->ConnectOutput(task, 0, mgr->CreateContainer(task->GetName(), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Task%s.root", task->GetName())));
+    if(TSTBIT(fSteerTask, kClusterErrorParam)){
+      mgr->AddTask(task = new AliTRDclusterResolution());
+      taskPtr[(Int_t)kClusterErrorParam] = task;
+      mgr->ConnectInput(task, 0, co[0]);
+      mgr->ConnectOutput(task, 0, mgr->CreateContainer(task->GetName(), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Task%s.root", task->GetName())));
+  
+      mgr->AddTask(task = new AliTRDclusterResolution());
+      mgr->ConnectInput(task, 0, co[1]);
+      mgr->ConnectOutput(task, 0, mgr->CreateContainer(Form("%sMC", task->GetName()), TObjArray::Class(), AliAnalysisManager::kOutputContainer, Form("TRD.Task%sMC.root", task->GetName())));
+    }
   }
 
   //____________________________________________
@@ -259,7 +271,7 @@ void run(Char_t *tasks="ALL", const Char_t *files=0x0, Int_t nmax=-1)
 
   if (!mgr->InitAnalysis()) return;
   printf("\n\tRUNNING TRAIN FOR TASKS:\n");
-  for(Int_t itask = 1; itask < NTRDTASKS; itask++){
+  for(Int_t itask = 1; itask <= NTRDTASKS; itask++){
     if(TSTBIT(fSteerTask, itask)) printf("\t   %s [%s]\n", taskPtr[itask]->GetTitle(), taskPtr[itask]->GetName());
   }
   printf("\n\n");

diff --git a/TRD/qaRec/run.h b/TRD/qaRec/run.h
index 1cc08e6c16d8f5e9ec53044955d86ca7bcc41171..87653fbe8e5147af3f2510bd15e81f00ae559381 100644 (file)
--- a/TRD/qaRec/run.h
+++ b/TRD/qaRec/run.h
@@ -6,7 +6,7 @@
 #define TSTBIT(n,i) ((Bool_t)(((n) & BIT(i)) != 0))
 #define CLRBIT(n,i) ((n) &= ~BIT(i))
 
-#define NTRDTASKS 7
+#define NTRDTASKS 8
 
 enum AliTRDrecoTasks{
    kInfoGen = 0
@@ -17,6 +17,7 @@ enum AliTRDrecoTasks{
   ,kCalibration = 5
   ,kPIDChecker = 6
   ,kPIDRefMaker = 7
+  ,kClusterErrorParam = 8
 };
 
 const Char_t* fgkTRDtaskClassName[NTRDTASKS] = {
@@ -27,6 +28,7 @@ const Char_t* fgkTRDtaskClassName[NTRDTASKS] = {
   ,"AliTRDcalibration"
   ,"AliTRDpidChecker"
   ,"AliTRDpidRefMaker"
+  ,"AliTRDclusterResolution"
 };
 
 const Char_t *fgkTRDtaskOpt[NTRDTASKS+3] = {
@@ -38,6 +40,7 @@ const Char_t *fgkTRDtaskOpt[NTRDTASKS+3] = {
   ,"CAL"
   ,"PID"
   ,"PIDR"
+  ,"CLRES"
   ,"NOFR"
   ,"NOMC"
 };