]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDQADataMakerRec.cxx
better ESD performance plot
[u/mrichter/AliRoot.git] / TRD / AliTRDQADataMakerRec.cxx
index eedcff18b22a3e0675a1a2cf65dc6a12da506588..855165e30b00ea604a9abdbb3a1c1a0c03969b18 100644 (file)
 ////////////////////////////////////////////////////////////////////////////
 
 // --- ROOT system ---
-#include <TClonesArray.h>
-#include <TFile.h> 
 #include <TH1D.h> 
 #include <TH2D.h>
 #include <TH3D.h>
 #include <TProfile.h>
 #include <TF1.h>
-#include <TCanvas.h>
 
 // --- AliRoot header files ---
 #include "AliESDEvent.h"
-#include "AliLog.h"
+#include "AliRawReader.h"
+#include "AliQAChecker.h"
+
 #include "AliTRDcluster.h"
 #include "AliTRDQADataMakerRec.h"
-#include "AliTRDgeometry.h"
-#include "AliTRDdataArrayI.h"
-#include "AliTRDrawStreamTB.h"
-
-#include "AliQAChecker.h"
+#include "AliTRDrawStream.h"
 
 ClassImp(AliTRDQADataMakerRec)
 
 //____________________________________________________________________________ 
   AliTRDQADataMakerRec::AliTRDQADataMakerRec() : 
-  AliQADataMakerRec(AliQA::GetDetName(AliQA::kTRD), "TRD Quality Assurance Data Maker")
+  AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kTRD), "TRD Quality Assurance Data Maker")
 {
   //
   // Default constructor
+  //
+
 }
 
 //____________________________________________________________________________ 
@@ -83,165 +80,175 @@ AliTRDQADataMakerRec& AliTRDQADataMakerRec::operator=(const AliTRDQADataMakerRec
 }
 
 //____________________________________________________________________________ 
-void AliTRDQADataMakerRec::EndOfDetectorCycle(AliQA::TASKINDEX task, TObjArray * list)
+void AliTRDQADataMakerRec::EndOfDetectorCycle(AliQAv1::TASKINDEX_t task, TObjArray ** list)
 {
   //
   // Detector specific actions at end of cycle
   //
   //TStopwatch watch;
   //watch.Start();
+  /**/
+  AliDebug(AliQAv1::GetQADebugLevel(), "End of TRD cycle");
   
-  AliInfo("End of TRD cycle");
-  
-  if (task == AliQA::kRECPOINTS) {
+  if (task == AliQAv1::kRECPOINTS) {
     
-    TH1D *hist = new TH1D("fitHist", "", 200, -0.5, 199.5);
-    //list->Print();
-    
-    // fill detector map;
-    for(int i=0; i<540; i++) {
-      Double_t v = ((TH1D*)list->At(0))->GetBinContent(i+1);
-      Int_t sm = i/30;
-      Int_t det = i%30;
-
-      TH2D *detMap = (TH2D*)list->At(87);
-      Int_t bin = detMap->FindBin(sm, det);
-      detMap->SetBinContent(bin, v);
-    }
+    TH1D * hist = new TH1D("fitHist", "", 200, -0.5, 199.5);
 
-
-    // Rec points full chambers
-    for (Int_t i=0; i<540; i++) {
+    // loop over event types
+    for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
+      if (! IsValidEventSpecie(specie, list)) 
+        continue ;
+    //list[specie]->Print();
       
-      //AliInfo(Form("I = %d", i));
-
-      //TH1D *h = ((TH2D*)list->At(1))->ProjectionY(Form("qaTRD_recPoints_amp_%d",i), i+1, i+1);
-      hist->Reset();
-      for(Int_t b=1; b<hist->GetXaxis()->GetNbins()-1; b++) {
-       Double_t xvalue = hist->GetBinCenter(b);
-       Int_t bin = ((TH2D*)list->At(1))->FindBin(i,xvalue);
-       Double_t value =  ((TH2D*)list->At(1))->GetBinContent(bin);
-       hist->SetBinContent(b, value);
+      // fill detector map;
+      for(Int_t i = 0 ; i < 540 ; i++) {
+        Double_t v = ((TH1D*)list[specie]->At(0))->GetBinContent(i+1);
+        Int_t sm = i/30;
+        Int_t det = i%30;
+        TH2D *detMap = (TH2D*)list[specie]->At(87);
+        Int_t bin = detMap->FindBin(sm, det);
+        detMap->SetBinContent(bin, v);
       }
       
-      //printf("Sum = %d %f\n", i, hist->GetSum());
-      if (hist->GetSum() < 100) continue; // chamber not present
-      
-      hist->Fit("landau", "q0", "goff", 10, 180);
-      TF1 *fit = hist->GetFunction("landau");
-      ((TH1D*)list->At(12))->Fill(fit->GetParameter(1));
-      ((TH1D*)list->At(13))->Fill(fit->GetParameter(2));
-    }
-    
-    // time-bin by time-bin sm by sm
-    for(Int_t i=0; i<18; i++) { // loop over super-modules
-      
-      for(Int_t j=0; j<35; j++) { // loop over time bins
-       
-       //TH1D *h =  ((TH3D*)list->At(10))->ProjectionZ(Form("ampTime_%d",i), i+1, i+1, j+1, j+1);     
-       hist->Reset();
-       for(Int_t b=1; b<hist->GetXaxis()->GetNbins()-1; b++) {
-         Double_t xvalue = hist->GetBinCenter(b);
-         Double_t svalue = 0;
-         
-         for(Int_t det=i*30; det<(i+1)*30; det++) { // loop over detectors
-           Int_t bin = ((TH3D*)list->At(10))->FindBin(det,j,xvalue);
-           Double_t value =  ((TH3D*)list->At(10))->GetBinContent(bin);
-           svalue += value;
-         }
-         //printf("v = %f\n", value);
-         hist->SetBinContent(b, svalue);
-       }
+      // Rec points full chambers
+      for (Int_t i = 0 ; i < 540 ; i++) {
+        //AliDebug(AliQAv1::GetQADebugLevel(), Form("I = %d", i));
+        //TH1D *h = ((TH2D*)list[specie]->At(1))->ProjectionY(Form("qaTRD_recPoints_amp_%d",i), i+1, i+1);
+        hist->Reset();
        
-       if (hist->GetSum() < 100) continue;
-       //printf("fitting %d %d %f\n", i, j, hist->GetSum());
+       // project TH2D into TH1D 
+        for(Int_t b = 1 ; b < hist->GetXaxis()->GetNbins()-1 ; b++) {
+          Double_t xvalue = hist->GetBinCenter(b);
+          Int_t bin = ((TH2D*)list[specie]->At(1))->FindBin(i,xvalue);
+          Double_t value =  ((TH2D*)list[specie]->At(1))->GetBinContent(bin);
+          hist->SetBinContent(b, value);
+        }
+        //AliDebug(AliQAv1::GetQADebugLevel(), Form("Sum = %d %f\n", i, hist->GetSum()));
+        if (hist->GetSum() < 100) 
+         continue; // not enougth data in a chamber
        
-       hist->Fit("landau", "q0", "goff", 10, 180);
-       TF1 *fit = hist->GetFunction("landau");
-       
-       TH1D *h1 = (TH1D*)list->At(14+18+i);
-       Int_t bin = h1->FindBin(j);
-       // printf("%d %d %d\n", det, j, bin);
-       h1->SetBinContent(bin, TMath::Abs(fit->GetParameter(1)));
+        hist->Fit("landau", "q0", "goff", 10, 180);
+        TF1 *fit = hist->GetFunction("landau");
+        ((TH1D*)list[specie]->At(12))->Fill(fit->GetParameter(1));
+        ((TH1D*)list[specie]->At(13))->Fill(fit->GetParameter(2));
       }
-    }
-      
 
-    // time-bin by time-bin chamber by chamber
-    for (Int_t i=0; i<540; i++) {
-      
-      //TH1D *test = ((TH3D*)list->At(10))->ProjectionZ(Form("ampTime_%d",i), i+1, i+1, 0, 35);     
-      //if (test->GetSum() < 100) continue;
-      
-      //AliInfo(Form("fitting det = %d", i));
-      
-      for(Int_t j=0; j<35; j++) {
-       
-       //TH1D *h =  ((TH3D*)list->At(10))->ProjectionZ(Form("ampTime_%d",i), i+1, i+1, j+1, j+1);     
-       hist->Reset();
-       for(Int_t b=1; b<hist->GetXaxis()->GetNbins()-1; b++) {
-         Double_t xvalue = hist->GetBinCenter(b);
-         Int_t bin = ((TH3D*)list->At(10))->FindBin(i,j,xvalue);
-         Double_t value =  ((TH3D*)list->At(10))->GetBinContent(bin);
-         //printf("v = %f\n", value);
-         hist->SetBinContent(b, value);
+
+      // time-bin by time-bin sm by sm
+      for(Int_t i=0; i<18; i++) { // loop over super-modules      
+       for(Int_t j=0; j<kTimeBin; j++) { // loop over time bins
+         
+          hist->Reset();
+          for(Int_t b = 1 ; b < hist->GetXaxis()->GetNbins()-1 ; b++) {
+            Double_t xvalue = hist->GetBinCenter(b);
+            Double_t svalue = 0.0;
+            for(Int_t det = i*30 ; det < (i+1)*30 ; det++) { // loop over detectors
+              Int_t bin = ((TH3D*)list[specie]->At(10))->FindBin(det,j,xvalue);
+              Double_t value =  ((TH3D*)list[specie]->At(10))->GetBinContent(bin);
+              svalue += value;
+            }
+            //AliDebug(AliQAv1::GetQADebugLevel(), Form("v = %f\n", value));
+            hist->SetBinContent(b, svalue);
+          }
+         
+          if (hist->GetSum() < 100) 
+            continue;
+         
+         hist->Fit("landau", "q0", "goff", 10, 180);
+         TF1 *fit = hist->GetFunction("landau");
+         
+         TH1D *h1 = (TH1D*)list[specie]->At(14+18+i);
+         h1->SetMarkerStyle(20);
+         Int_t bin = h1->FindBin(j);
+         // printf("%d %d %d\n", det, j, bin);
+         
+         Double_t value = TMath::Abs(fit->GetParameter(1));
+         Double_t error = TMath::Abs(fit->GetParError(1));
+        
+         if (value/error < 3) continue; // insuficient statistics
+         
+         h1->SetBinContent(bin, value);
+         h1->SetBinError(bin, error);  
        }
+      }
        
-       if (hist->GetSum() < 100) continue;
-       //printf("fitting %d %d %f\n", i, j, hist->GetSum());
-       
-       hist->Fit("landau", "q0", "goff", 10, 180);
-       TF1 *fit = hist->GetFunction("landau");
+      // for numerical convergence
+      TF1 *form = new TF1("formLandau", "landau", 0, 200);
        
-       Int_t sm = i/30;
-       Int_t det = i%30;
-       TH2D *h2 = (TH2D*)list->At(14+sm);
-       Int_t bin = h2->FindBin(det,j);
-       // printf("%d %d %d\n", det, j, bin);
-       h2->SetBinContent(bin, TMath::Abs(fit->GetParameter(1)));
-       h2->SetBinError(bin,fit->GetParError(1));
+      // time-bin by time-bin chamber by chamber
+      for (Int_t i=0; i<540; i++) {
+       for(Int_t j=0; j<kTimeBin; j++) {
+           
+         hist->Reset();
+         for(Int_t b = 1 ; b < hist->GetXaxis()->GetNbins()-1 ; b++) {
+           Double_t xvalue = hist->GetBinCenter(b);
+           Int_t bin = ((TH3D*)list[specie]->At(10))->FindBin(i,j,xvalue);
+           Double_t value =  ((TH3D*)list[specie]->At(10))->GetBinContent(bin);
+           //AliDebug(AliQAv1::GetQADebugLevel(), Form("v = %f\n", value));
+           hist->SetBinContent(b, value);
+         }
+         
+         if (hist->GetSum() < 100) 
+           continue;
+         
+         form->SetParameters(1000, 60, 20);
+         hist->Fit(form, "q0", "goff", 20, 180);
+         
+         Int_t sm = i/30;
+         Int_t det = i%30;
+         TH2D *h2 = (TH2D*)list[specie]->At(14+sm);
+         Int_t bin = h2->FindBin(det,j);
+         // printf("%d %d %d\n", det, j, bin);
+         
+         Double_t value = TMath::Abs(form->GetParameter(1));
+         Double_t error = TMath::Abs(form->GetParError(1));
+         
+         if (value/error < 3) continue;
+         
+         h2->SetBinContent(bin, value);
+         h2->SetBinError(bin, error);
+       }
       }
     }
-    
-    if (hist) delete hist;
+    if (hist) 
+      delete hist;
   }
-  
   //////////////////////////
   // const Int_t knbits = 6;
   // const char *suf[knbits] = {"TPCi", "TPCo", "TPCz", "TRDo", "TRDr", "TRDz"};
   //const char *sufRatio[4] = {"TRDrTRDo", "TRDoTPCo", "TRDrTPCo", "TRDzTPCo"};
 
-  if (task == AliQA::kESDS) {
+  if (task == AliQAv1::kESDS) {
     
     const Int_t knRatio = 4;
     const Int_t kN[knRatio] = {4,3,4,5};
     const Int_t kD[knRatio] = {3,1,1,3}; 
     
     // create ratios
-    for(Int_t type=0; type<2; type++) {
-      for(Int_t i=0; i<knRatio; i++) {
-
-       TH1D *ratio = (TH1D*)list->At(19 + 2*i + type);
-       TH1D *histN = (TH1D*)list->At(3 + 2*kN[i] + type);
-       TH1D *histD = (TH1D*)list->At(3 + 2*kD[i] + type);
-
-       BuildRatio(ratio, histN, histD);
-       //ratio->Reset();
-       //ratio->Add(histN);
-       //ratio->Divide(histD);
+    for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
+      if (! IsValidEventSpecie(specie, list)) 
+        continue ;
+     for(Int_t type = 0 ; type < 2 ; type++) {
+        for(Int_t i = 0 ; i < knRatio ; i++) {
+          TH1D *ratio = (TH1D*)list[specie]->At(19 + 2*i + type);
+          TH1D *histN = (TH1D*)list[specie]->At(3 + 2*kN[i] + type);
+          TH1D *histD = (TH1D*)list[specie]->At(3 + 2*kD[i] + type);
+          BuildRatio(ratio, histN, histD);
+          //ratio->Reset();
+          //ratio->Add(histN);
+          //ratio->Divide(histD);
+        }
       }
+      // ratio for the fraction of electrons per stack
+      TH1D *histN = (TH1D*)list[specie]->At(33);
+      TH1D *histD = (TH1D*)list[specie]->At(32);
+      TH1D *ratio = (TH1D*)list[specie]->At(34);
+      BuildRatio(ratio, histN, histD);
     }
-
-    // ratio for the fraction of electrons per stack
-    TH1D *histN = (TH1D*)list->At(33);
-    TH1D *histD = (TH1D*)list->At(32);
-    TH1D *ratio = (TH1D*)list->At(34);
-    BuildRatio(ratio, histN, histD);
   }
-
-  
   // call the checker
-  AliQAChecker::Instance()->Run(AliQA::kTRD, task, list) ;    
+  AliQAChecker::Instance()->Run(AliQAv1::kTRD, task, list) ;    
+
 }
 
 //____________________________________________________________________________ 
@@ -250,31 +257,33 @@ void AliTRDQADataMakerRec::InitESDs()
   //
   // Create ESDs histograms in ESDs subdir
   //
-
+  const Bool_t expert   = kTRUE ; 
+  const Bool_t image    = kTRUE ; 
+  
   const Int_t kNhist = 36+5+4;
 
   TH1 *hist[kNhist];
   Int_t histoCounter = -1 ;
 
-  hist[++histoCounter] = new TH1D("qaTRD_esd_ntracks", ";Number of tracks", 300, -0.5, 299.5);
-  hist[++histoCounter] = new TH1D("qaTRD_esd_sector", ";Sector", 18, -0.5, 17.7);
-  hist[++histoCounter] = new TH1D("qaTRD_esd_bits", ";Bits", 64, -0.5, 63.5);
+  hist[++histoCounter] = new TH1D("qaTRD_esd_ntracks", "TRD esd ntracks;Number of tracks;Counts", 300, -0.5, 299.5);
+  hist[++histoCounter] = new TH1D("qaTRD_esd_sector", "TRD esd sector;Sector;Counts", 18, -0.5, 17.7);
+  hist[++histoCounter] = new TH1D("qaTRD_esd_bits", "TRD esd bits;Bits;Counts", 64, -0.5, 63.5);
 
   const Int_t knbits = 6;
   const char *suf[knbits] = {"TPCi", "TPCo", "TPCz", "TRDo", "TRDr", "TRDz"};
 
   // histo = 3
   for(Int_t i=0; i<knbits; i++) { 
-    hist[++histoCounter] = new TH1D(Form("qaTRD_esd_pt%s",suf[i]), ";p_{T} (GeV/c);", 100, 0, 10);
+    hist[++histoCounter] = new TH1D(Form("qaTRD_esd_pt%s",suf[i]), Form("qaTRD_esd_pt%s;p_{T} (GeV/c);Counts",suf[i]), 100, 0, 10);
     hist[++histoCounter] = new TH1D(Form("qaTRD_esd_trdz%s", suf[i]), ";z (cm)", 200, -400, 400); 
   }
 
-  hist[++histoCounter] = new TH1D("qaTRD_esd_clsTRDo", "TRDo;number of clusters", 180, -0.5, 179.5);;
-  hist[++histoCounter] = new TH1D("qaTRD_esd_clsTRDr", "TRDr;number of clusters", 180, -0.5, 179.5);;
-  hist[++histoCounter] = new TH1D("qaTRD_esd_clsTRDz", "TRDz;number of clusters", 180, -0.5, 179.5);;
+  hist[++histoCounter] = new TH1D("qaTRD_esd_clsTRDo", "TRDo;number of clusters;Counts", 180, -0.5, 179.5);;
+  hist[++histoCounter] = new TH1D("qaTRD_esd_clsTRDr", "TRDr;number of clusters;Counts", 180, -0.5, 179.5);;
+  hist[++histoCounter] = new TH1D("qaTRD_esd_clsTRDz", "TRDz;number of clusters;Counts", 180, -0.5, 179.5);;
   //hist[++histoCounter] = new TH1D("qaTRD_esd_clsRatio", ";cluster ratio", 100, 0., 1.3);;
 
-  hist[++histoCounter] = new TH2D("qaTRD_esd_sigMom", ";momentum (GeV/c);signal", 100, 0, 5, 200, 0, 1e3);
+  hist[++histoCounter] = new TH2D("qaTRD_esd_sigMom", "TRD esd sig Mom;momentum (GeV/c);signal", 100, 0, 5, 200, 0, 1e3);
 
   // end of cycle plots (hist 19)
   const char *sufRatio[4] = {"TRDrTRDo", "TRDoTPCo", "TRDrTPCo", "TRDzTPCo"};
@@ -290,49 +299,51 @@ void AliTRDQADataMakerRec::InitESDs()
   }
 
   // 27 - 31
-  hist[27] = new TH1D("qaTRD_esd_quality", ";quality", 120, 0, 12);
-  hist[28] = new TH1D("qaTRD_esd_budget", ";NN", 110, -1000, 100);
-  hist[29] = new TH1D("qaTRD_esd_chi2", ";chi2", 200, 0, 100);
-  hist[30] = new TH1D("qaTRD_esd_timeBin", ";time bin", 7, -0.5, 6.5);
-  hist[31] = new TH1D("qaTRD_esd_pidQuality", "pid Quality", 7, -0.5, 6.5);
+  hist[27] = new TH1D("qaTRD_esd_quality", "TRD esd quality;quality;Counts", 120, 0, 12);
+  hist[28] = new TH1D("qaTRD_esd_budget", "TRD esd budget;NN;Counts", 110, -1000, 100);
+  hist[29] = new TH1D("qaTRD_esd_chi2", "TRD esd chi2;chi2;Counts", 200, 0, 100);
+  hist[30] = new TH1D("qaTRD_esd_timeBin", "TRD esd timeBin;time bin;Counts", 7, -0.5, 6.5);
+  hist[31] = new TH1D("qaTRD_esd_pidQuality", "pid Quality;quality;;Counts", 7, -0.5, 6.5);
 
   // stack by stack electron identyfication
-  hist[32] = new TH1D("qaTRD_esd_tracksStack", "number of all tracks;stack", 90, -0.5, 89.5);
-  hist[33] = new TH1D("qaTRD_esd_electronStack", "number of electron tracks;stack", 90, -0.5, 89.5);
-  hist[34] = new TH1D("qaTRD_esd_elRatioStack", "fraction of electron tracks;stack", 90, -0.5, 89.5);
+  hist[32] = new TH1D("qaTRD_esd_tracksStack", "number of all tracks;stack;Counts", 90, -0.5, 89.5);
+  hist[33] = new TH1D("qaTRD_esd_electronStack", "number of electron tracks;stack;Counts", 90, -0.5, 89.5);
+  hist[34] = new TH1D("qaTRD_esd_elRatioStack", "fraction of electron tracks;stack;Counts", 90, -0.5, 89.5);
   hist[35] = new TH1D("qaTRD_esd_thetaOut", ";tan(theta);", 100, -1, 1);
   
   const char *partType[5] = {"Electron", "Muon", "Pion", "Kaon", "Proton"}; 
 
   for(Int_t i=0; i<AliPID::kSPECIES; i++)
     hist[36+i] = new TH1D(Form("qaTRD_esd_pid%d",i),
-                         Form("%s;probability",partType[i]), 100, 0, 1);
+                         Form("%s;probability;Counts",partType[i]), 100, 0, 1);
  
   // dE/dX vs momentum in three regions
   const char *zoneName[4] = {"total charge", "ampilification range", "plateau", "TR range"};
  
   // prepare the scale from 0.1 to 10 GeV
-  Double_t scalex[101];
-  Double_t dd = (TMath::Log(10) - TMath::Log(0.5)) / 100.;
-  for(Int_t ix=0; ix<101; ix++) {
+  const Int_t nscalex= 50;
+  Double_t scalex[nscalex+1];
+  Double_t dd = (TMath::Log(10) - TMath::Log(0.5)) / nscalex;
+  for(Int_t ix=0; ix<nscalex+1; ix++) {
     scalex[ix] = 0.5 * TMath::Exp(dd * ix);
   }
 
-  Double_t scaley[101];
-  for(Int_t iy=0; iy<101; iy++) {
-    scaley[iy] = iy * (3e3/100.);
+  const Int_t nscaley = 50;
+  Double_t scaley[nscaley+1];
+  for(Int_t iy=0; iy<nscaley+1; iy++) {
+    scaley[iy] = iy * (3e3/nscaley);
   }
     
 
   for(Int_t i=0; i<4; i++) {
     hist[41+i] = new TH2D(Form("qaTRD_esd_signalPzone_%d",i), 
-                         Form("%s;momentum (GeV/c);singal (a.u.)", zoneName[i]),
-                         100, scalex, 100, scaley);
+                         Form("%s;momentum (GeV/c);signal (a.u.)", zoneName[i]),
+                         nscalex, scalex, nscaley, scaley);
   }
 
   for(Int_t i=0; i<kNhist; i++) {
     //hist[i]->Sumw2();
-    Add2ESDsList(hist[i], i);
+    Add2ESDsList(hist[i], i, !expert, image);
   }
 
 }
@@ -343,69 +354,84 @@ void AliTRDQADataMakerRec::InitRecPoints()
   //
   // Create Reconstructed Points histograms in RecPoints subdir
   //
+  const Bool_t expert   = kTRUE ; 
+  const Bool_t image    = kTRUE ; 
+  
+  //printf("Helo from Init rec points\n");
 
-  const Int_t kNhist = 14 + 4 * 18 + 2;
+  const Int_t kNhist = 14 + 4 * 18 + 2 + 9;// + 540;
   TH1 *hist[kNhist];
 
-  hist[0] = new TH1D("qaTRD_recPoints_det", ";Detector ID of the cluster", 540, -0.5, 539.5);
-  hist[1] = new TH2D("qaTRD_recPoints_amp", ";Amplitude", 540, -0.5, 539, 200, -0.5, 199.5);
-  hist[2] = new TH1D("qaTRD_recPoints_npad", ";Number of Pads", 12, -0.5, 11.5);
+  hist[0] = new TH1D("qaTRD_recPoints_det", "RRD recPoints det;Detector ID of the cluster;Counts", 540, -0.5, 539.5);
+  hist[1] = new TH2D("qaTRD_recPoints_amp", "TRD recPoints amp;Amplitude;??", 540, -0.5, 539, 200, -0.5, 199.5);
+  hist[2] = new TH1D("qaTRD_recPoints_npad", "TRD recPoints npad;Number of Pads;Counts", 12, -0.5, 11.5);
 
-  hist[3] = new TH1D("qaTRD_recPoints_dist2", ";residuals [2pad]", 100, -1, 1);
-  hist[4] = new TH1D("qaTRD_recPoints_dist3", ";residuals [3pad]", 100, -1, 1);
-  hist[5] = new TH1D("qaTRD_recPoints_dist4", ";residuals [4pad]", 100, -1, 1);
-  hist[6] = new TH1D("qaTRD_recPoints_dist5", ";residuals [5pad]", 100, -1, 1);
+  hist[3] = new TH1D("qaTRD_recPoints_dist2", "TRD recPoints dist2;residuals [2pad];Counts", 100, -1, 1);
+  hist[4] = new TH1D("qaTRD_recPoints_dist3", "TRD recPoints dist3;residuals [3pad];Counts", 100, -1, 1);
+  hist[5] = new TH1D("qaTRD_recPoints_dist4", "TRD recPoints dist4;residuals [4pad];Counts", 100, -1, 1);
+  hist[6] = new TH1D("qaTRD_recPoints_dist5", "TRD recPoints dist5;residuals [5pad];Counts", 100, -1, 1);
 
-  hist[7] = new TH2D("qaTRD_recPoints_rowCol", ";row;col", 16, -0.5, 15.5, 145, -0.5, 144.5);
-  hist[8] = new TH1D("qaTRD_recPoints_time", ";time bin", 35, -0.5, 34.5);
-  hist[9] = new TH1D("qaTRD_recPoints_nCls", ";number of clusters", 500, -0.5, 499.5);
+  hist[7] = new TH2D("qaTRD_recPoints_rowCol", "TRDrecPointsrowCol;row;col", 16, -0.5, 15.5, 145, -0.5, 144.5);
+  hist[8] = new TH1D("qaTRD_recPoints_time", "TRDrecPoints time;time bin;Counts", kTimeBin, -0.5, kTimeBin-0.5);
+  hist[9] = new TH1D("qaTRD_recPoints_nCls", "TRD recPoints nCls;number of clusters;Counts", 500, -0.5, 499.5);
 
-  hist[10] = new TH3D("qaTRD_recPoints_sigTime", ";chamber;time bin;signal", 
-                     540, -0.5, 539.5, 35, -0.5, 34.5, 200, -0.5, 199.5);
-  hist[11] = new TProfile("qaTRD_recPoints_prf", ";distance;center of gravity"
+  hist[10] = new TH3D("qaTRD_recPoints_sigTime", "TRD recPoints sigTime;chamber;time bin;signal", 
+                     540, -0.5, 539.5, kTimeBin, -0.5, kTimeBin-0.5, 200, -0.5, 199.5);
+  hist[11] = new TProfile("qaTRD_recPoints_prf", "TRD recPoints prf;distance;center of gravity;Counts"
                          , 120, -0.6, 0.6, -1.2, 1.2, "");
 
-  hist[12] = new TH1D("qaTRD_recPoints_ampMPV", ";amplitude MPV", 200, 0, 200);
-  hist[13] = new TH1D("qaTRD_recPoints_ampSigma", ";amplitude Sigma", 200, 0, 200); 
+  hist[12] = new TH1D("qaTRD_recPoints_ampMPV", "TRD recPoints ampMPV;amplitude MPV;Counts", 150, 0, 150);
+  hist[13] = new TH1D("qaTRD_recPoints_ampSigma", "TRD recPoints ampSigma;amplitude Sigma;Counts", 200, 0, 200); 
   
   // chamber by chamber
   for(Int_t i=0; i<18; i++) {
     hist[14+i] = new TH2D(Form("qaTRD_recPoints_sigTime_sm%d",i), Form("sm%d;det;time bin"), 
-                       30, -0.5, 29.5, 35, -0.5, 34.5);
-    hist[14+i]->SetMinimum(20);
-    hist[14+i]->SetMaximum(40);
+                       30, -0.5, 29.5, kTimeBin, -0.5, kTimeBin-0.5);
+    hist[14+i]->SetMinimum(0);
+    hist[14+i]->SetMaximum(150);
   }
  
   // time bin by time bin sm-by-sm
   for(Int_t i=0; i<18; i++) {
     hist[14+18+i] = new TH1D(Form("qaTRD_recPoints_sigTimeShape_sm%d", i), 
                             Form("sm%d;time bin;signal"),
-                            35, -0.5, 34.5);
+                            kTimeBin, -0.5, kTimeBin-0.5);
 
-    hist[14+18+i]->SetMaximum(120);    
+    hist[14+18+i]->SetMaximum(150);    
   }
 
   // str = 50
   for(Int_t i=0; i<18; i++) {
     hist[50+i] = new TH1D(Form("qaTRD_recPoints_nCls_sm%d",i),
                          Form("sm%d;time bin;number of clusters",i),
-                         35, -0.5, 34.5);
+                         kTimeBin, -0.5, kTimeBin-0.5);
   }
 
   // str = 68
   for(Int_t i=0; i<18; i++) {
     hist[68+i] = new TH1D(Form("qaTRD_recPoints_totalCharge_sm%d", i),
                          Form("sm%d;time bin;total charge", i),
-                         35, -0.5, 34.5);
+                         kTimeBin, -0.5, kTimeBin-0.5);
   }
 
-  hist[86] = new TH1D("qaTRD_recPoints_signal", ";amplitude", 200, -0.5, 199.5);
-  hist[87] = new TH2D("qaTRD_recPoints_detMap", ";sm;chamber", 18, -0.5, 17.5, 30, -0.5, 29.5);
+  hist[86] = new TH1D("qaTRD_recPoints_signal", "TRD recPoints signal;amplitude;Counts", 400, -0.5, 399.5);
+  hist[87] = new TH2D("qaTRD_recPoints_detMap", "TRD recPoints detMap;sm;chamber;Counts", 18, -0.5, 17.5, 30, -0.5, 29.5);
+
+  
+  // amplitude as a function of the pad size
+  for(Int_t i=0; i<9; i++) {
+    hist[88+i] = new TH1D(Form("qaTRD_recPoints_signalNpad_%d", i+2), Form("qaTRD_recPoints_signalNpad_%d;amplitude, ADC", i+2), 400, -0.5, 399.5); 
+  }
+  
+  // one 2D histogram per chamber
+  //  for(Int_t i=0; i<540; i++) {
+  //  hist[88+i] = new TH2D(Form("qaTRD_recPoints_map%d", i), ";col;row", 16, -0.5, 15.5, 144, -0.5, 143.5);
+  //}
 
 
   for(Int_t i=0; i<kNhist; i++) {
     //hist[i]->Sumw2();
-    Add2RecPointsList(hist[i], i);
+    Add2RecPointsList(hist[i], i, !expert, image);
   }
 }
 
@@ -415,12 +441,30 @@ void AliTRDQADataMakerRec::InitRaws()
   //
   // create Raws histograms in Raws subdir
   //
+  const Bool_t expert   = kTRUE ; 
+  const Bool_t saveCorr = kTRUE ; 
+  const Bool_t image    = kTRUE ; 
+  
+  AliInfo("Initialization of QA for Raw Data");
 
-  const Int_t kSM = 18;
+  //const Int_t kSM = 18;
   //const Int_t kNCh = 540;
-  const Int_t kNhist = 4+kSM;
-  TH1D *hist[kNhist];
+  const Int_t kNhist = 6;
+  TH1 *hist[kNhist];
+
+  // link monitor
+  const char *linkName[3] = {"smLink", "smBeaf", "smData"};
+  for(Int_t i=0; i<3; i++) {
+    hist[i] = new TH2D(Form("qaTRD_raws_%s", linkName[i]), 
+                      ";super module;link", 
+                      18, -0.5, 17.5, 60, -0.5, 59.5);
+  }
+
+  hist[3] = new TH1D("qaTRD_raws_errorHC", "TRD raws error HC;error ID;Counts", 18, -3.5, 14.5);
+  hist[4] = new TH1D("qaTRD_raws_errorMCM", "TRD raws error MCM;error ID;Counts", 18, -3.5, 14.5);
+  hist[5] = new TH1D("qaTRD_raws_errorADC", "TRD raws errorADC;error ID;Counts", 18, -3.5, 14.5);
 
+  /*
   // four histograms to be published
   hist[0] = new TH1D("qaTRD_raws_det", ";detector", 540, -0.5, 539.5);
   hist[1] = new TH1D("qaTRD_raws_sig", ";signal", 100, -0.5, 99.5);
@@ -432,29 +476,30 @@ void AliTRDQADataMakerRec::InitRaws()
   const Int_t kNMCM = 30 * 8 * 16;
   for(Int_t i=0; i<kSM; i++)
     hist[4+i] = new TH1D(Form("qaTRD_raws_sm%d",i),"",kNMCM, -0.5, kNMCM-0.5); 
-  
+  */
+
   // register
   for(Int_t i=0; i<kNhist; i++) {
     //hist[i]->Sumw2();
-    Add2RawsList(hist[i], i);
+    Add2RawsList(hist[i], i, !expert, image, !saveCorr);
   }
 
 }
 
 //____________________________________________________________________________
-void AliTRDQADataMakerRec::MakeESDs(AliESDEvent * esd)
+void AliTRDQADataMakerRec::MakeESDs(AliESDEvent * const esd)
 {
   //
   // Make QA data from ESDs
   //
-
+  
   Int_t nTracks = esd->GetNumberOfTracks();
   GetESDsData(0)->Fill(nTracks);
 
   // track loop
-  for (Int_t i=0; i<nTracks; i++) {
+  for (Int_t iTrack = 0; iTrack<nTracks; iTrack++) {
 
-    AliESDtrack *track = esd->GetTrack(i);
+    AliESDtrack *track = esd->GetTrack(iTrack);
     const AliExternalTrackParam *paramOut = track->GetOuterParam();
     const AliExternalTrackParam *paramIn = track->GetInnerParam();
 
@@ -515,26 +560,26 @@ void AliTRDQADataMakerRec::MakeESDs(AliESDEvent * esd)
     GetESDsData(28)->Fill(track->GetTRDBudget());
     GetESDsData(29)->Fill(track->GetTRDchi2());
     GetESDsData(30)->Fill(track->GetTRDTimBin(0));
-    GetESDsData(31)->Fill(track->GetTRDpidQuality());
+    GetESDsData(31)->Fill(track->GetTRDntrackletsPID());
     
     
     // dedx
-    for(Int_t i=0; i<4; i++) {
+    for(Int_t k=0; k<4; ++k) {
       Double_t dedx = 0;
       for(Int_t j=0; j<6; j++) {
-       dedx += track->GetTRDsignals(j, i-1);
+       dedx += track->GetTRDslice(j, k-1);
       }
-      GetESDsData(41+i)->Fill(paramOut->GetP(), dedx/6.);
+      GetESDsData(41+k)->Fill(paramOut->GetP(), dedx/6.);
     }
 
     // probabilities
     if (status & AliESDtrack::kTRDpid) {
-      for(Int_t i=0; i<AliPID::kSPECIES; i++
-       GetESDsData(36+i)->Fill(track->GetTRDpid(i));
+      for(Int_t k=0; k<AliPID::kSPECIES; ++k
+       GetESDsData(36+k)->Fill(track->GetTRDpid(k));
     }
 
     // probabilities uniformity
-    if (track->GetTRDpidQuality() < 6) continue;
+    if (track->GetTRDntrackletsPID() < 6) continue;
     GetESDsData(35)->Fill(paramOut->GetZ()/paramOut->GetX());
     
     Int_t idx = 5 * sector + stack;
@@ -611,11 +656,11 @@ Int_t AliTRDQADataMakerRec::GetStack(const AliExternalTrackParam *paramOut) cons
   // calculates the stack the track is in
   //
   
-  const Double_t L = -0.9;
-  const Double_t W = (2*L)/5;
+  const Double_t l = -0.9;
+  const Double_t w = (2*l)/5;
 
   Double_t tan = paramOut->GetZ() / paramOut->GetX();
-  Double_t pos = (tan - L) / W;
+  Double_t pos = (tan - l) / w;
   return (Int_t) pos;
 }
 
@@ -645,7 +690,121 @@ void AliTRDQADataMakerRec::MakeRaws(AliRawReader* rawReader)
   //
   // Makes QA data from raw data
   //
+  
+  AliInfo("Execution of QA for Raw data");
+  //printf("Execution of QA for Raw data\n");
+
+  // create raw reader TB
+  rawReader->Reset();
+  rawReader->SelectEquipment(0, 1024, 1041);
+  rawReader->Select("TRD");
+  
+  // AliTRDrawStreamTB::AllowCorruptedData();
+  // AliTRDrawStreamTB::DisableStackNumberChecker();
+  // AliTRDrawStreamTB::DisableStackLinkNumberChecker();
+  // AliTRDrawStreamTB::DisableSkipData();
 
+  //AliTRDrawStreamTB *data = (AliTRDrawStreamTB*)AliTRDrawStreamBase::GetRawStream(rawReader);
+  AliTRDrawStreamBase::SetRawStreamVersion ("REAL") ;
+  AliTRDrawStream *data = (AliTRDrawStream *)AliTRDrawStreamBase::GetRawStream(rawReader);  
+
+  //if (raw->IsA()->GetName())
+
+
+  // import the histograms
+  TH2D *fSMLink[3];
+  for(Int_t i=0; i<3; i++) {
+    fSMLink[i] = (TH2D*)GetRawsData(i);
+    //printf("address = %d\n", fSMLink[i]);
+  }
+    
+  TH1D *fErrorHC =  (TH1D*)GetRawsData(3);
+  TH1D *fErrorMCM = (TH1D*)GetRawsData(4);
+  TH1D *fErrorADC = (TH1D*)GetRawsData(5);
+
+
+  // loop over super-modules
+  while (data->NextBuffer()>0) {
+    
+    //printf("processing next buffer\n");
+    // check sm
+    Int_t sm = rawReader->GetEquipmentId() - 1024;
+    if (sm < 0 || sm > 18) return;
+    
+    // loop over links
+    for (Int_t istack = 0; istack < 5; istack++) {     
+      for (Int_t ilink = 0; ilink < 12; ilink++) {
+       
+       //Int_t det = sm * 30 + istack * 6 + ilink/2;   
+      
+       // check if data delivered
+       if (!(data->IsLinkActiveInStack(istack, ilink))) continue;
+       fSMLink[0]->Fill(sm, istack * 12 + ilink);
+       
+       // check if beaf-beaf
+       if (data->GetLinkMonitorError(istack, ilink)) {
+         fSMLink[1]->Fill(sm, istack * 12 + ilink);
+         continue;
+       }
+       
+       // fill histogram with HC header errors
+       Int_t nErrHc = 0;
+       
+       nErrHc = FillBits(fErrorHC, data->GetH0ErrorCode(istack, ilink), 0);
+       if (!nErrHc) fErrorHC->Fill(-3);
+             
+       nErrHc = FillBits(fErrorHC, data->GetH1ErrorCode(istack, ilink), 2);
+       if (!nErrHc) fErrorHC->Fill(-2);
+             
+       nErrHc = FillBits(fErrorHC, data->GetHCErrorCode(istack, ilink), 4);
+       if (!nErrHc) fErrorHC->Fill(-1);
+       
+       // data integrity protection 
+       if (data->GetH0ErrorCode(istack, ilink) > 0) continue;
+       if (data->GetH1ErrorCode(istack, ilink) > 0) continue;
+       
+       fSMLink[2]->Fill(sm, istack * 12 + ilink);      
+
+       // loop over MCMs
+       for (Int_t imcm = 0; imcm < data->GetHCMCMmax(istack, ilink); imcm++ ){
+         
+         Int_t nErrMcm = 0;
+       
+         nErrMcm = FillBits(fErrorMCM, data->GetMCMhdErrorCode(istack, ilink, imcm), 0);
+         if (!nErrMcm) fErrorMCM->Fill(-3);
+       
+         nErrMcm = FillBits(fErrorMCM, data->GetMCMADCMaskErrorCode(istack, ilink, imcm), 5);
+         if (!nErrMcm) fErrorMCM->Fill(-2);
+       
+         nErrMcm = FillBits(fErrorMCM, data->GetMCMErrorCode(istack, ilink, imcm), 10);
+         if (!nErrMcm) fErrorMCM->Fill(-1);
+       
+         // MCM protection
+         if ( (data->GetMCMhdErrorCode(istack,ilink,imcm)) & 2 ) continue;
+
+         // loop over ADC chanels      
+         for (Int_t iadc=0; iadc < data->GetADCcount(istack, ilink, imcm); iadc++) {
+         
+           // fill ADC error bits
+           Int_t nErrAdc = FillBits(fErrorADC, data->GetADCErrorCode(), 0);
+           if (!nErrAdc) fErrorADC->Fill(-1);
+         }
+       } // mcm 
+      }
+    } // link
+
+    //printf("buffer analyzed\n");
+  }
+  
+  // clean up
+
+  delete data;
+  AliInfo("sucessfull execution of QA for TRD raw data");
+  //printf("sucessfull execution of QA for TRD raw data\n");
+  
+
+  /*
   // 157
   // T9 -- T10
 
@@ -657,9 +816,10 @@ void AliTRDQADataMakerRec::MakeRaws(AliRawReader* rawReader)
   const Int_t kMCM = 16;
   //  const Int_t kADC = 22;
 
+  rawReader->Reset() ; 
   //AliTRDrawStreamBase::SetRawStreamVersion("TB");
   AliTRDrawStreamBase *raw = AliTRDrawStreamBase::GetRawStream(rawReader);
-  AliInfo(Form("Stream version: %s", raw->IsA()->GetName()));
+  AliDebug(2,Form("Stream version: %s", raw->IsA()->GetName()));
 
   while (raw->Next()) {
 
@@ -686,6 +846,7 @@ void AliTRDQADataMakerRec::MakeRaws(AliRawReader* rawReader)
   }
 
   delete raw;
+  */
 }
 
 //____________________________________________________________________________
@@ -696,7 +857,7 @@ void AliTRDQADataMakerRec::MakeRecPoints(TTree * clustersTree)
   // 
   
   //  Info("MakeRecPoints", "making");
-
   Int_t nsize = Int_t(clustersTree->GetTotBytes() / (sizeof(AliTRDcluster))); 
   TObjArray *clusterArray = new TObjArray(nsize+1000); 
 
@@ -708,32 +869,68 @@ void AliTRDQADataMakerRec::MakeRecPoints(TTree * clustersTree)
   branch->SetAddress(&clusterArray); 
 
   // Loop through all entries in the tree
-  Int_t nEntries   = (Int_t) clustersTree->GetEntries();
+  Int_t nEntries   = (Int_t)TMath::Ceil( clustersTree->GetEntries() );
   Int_t nbytes     = 0;
   AliTRDcluster *c = 0;
   Int_t nDet[540];
   for (Int_t i=0; i<540; i++) nDet[i] = 0;
 
+  Int_t nCls = 0;
+  
+  //printf("nEntries = %d\n", nEntries);
+  //nEntries++;
+  
+  /*
+  // select the event 
+  for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
+
+    // Import the tree
+    nbytes += clustersTree->GetEvent(iEntry);  
+    Int_t nCluster = clusterArray->GetEntries();  
+    
+    for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
+      c = (AliTRDcluster *) clusterArray->At(iCluster);
+      nCls++;
+    }
+  }
+
+  if (nCls < 100) {
+    delete clusterArray;
+    return;
+  }
+  */
+
+  /////
+
   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
 
     // Import the tree
     nbytes += clustersTree->GetEvent(iEntry);  
 
     // Get the number of points in the detector
-    Int_t nCluster = clusterArray->GetEntriesFast();  
+    Int_t nCluster = clusterArray->GetEntries();  
 
+  
     // Loop through all TRD digits
     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
-      c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
+
+      nCls++;
+      c = (AliTRDcluster *) clusterArray->At(iCluster);
 
       Int_t iDet = c->GetDetector();
+      Int_t nPads = c->GetNPads();
+
       nDet[iDet]++;
       GetRecPointsData(0)->Fill(iDet);
       GetRecPointsData(86)->Fill(c->GetQ());
       GetRecPointsData(1)->Fill(iDet, c->GetQ());
-      GetRecPointsData(2)->Fill(c->GetNPads());
-      if (c->GetNPads() < 6)
+      GetRecPointsData(2)->Fill(nPads);
+      if (nPads < 6)
        GetRecPointsData(1+c->GetNPads())->Fill(c->GetCenter());
+      
+      if (nPads < 10)
+       GetRecPointsData(88+nPads-2)->Fill(c->GetQ());
+      else GetRecPointsData(96)->Fill(c->GetQ());
 
       //if (c->GetPadTime() < 5)
       ((TH2D*)GetRecPointsData(7))->Fill(c->GetPadRow(), c->GetPadCol());
@@ -744,6 +941,10 @@ void AliTRDQADataMakerRec::MakeRecPoints(TTree * clustersTree)
       Int_t iSM = iDet / 30;
       GetRecPointsData(50+iSM)->Fill(c->GetPadTime());
       GetRecPointsData(68+iSM)->Fill(c->GetPadTime(), c->GetQ());
+      
+      // total charge sm / det / timeBin
+      //((TH2D*)GetRecPointsData(14+iSM))->Fill(iDet-iSM*30, c->GetPadTime(), c->GetQ());
+
 
       // PRF for 2pad
       //if (c->GetNPads() == 2) {
@@ -762,9 +963,12 @@ void AliTRDQADataMakerRec::MakeRecPoints(TTree * clustersTree)
     }
   }
 
+  /*
   for(Int_t i=0; i<540; i++) 
     if (nDet[i] > 0) GetRecPointsData(9)->Fill(nDet[i]);
-
+  */
+  GetRecPointsData(9)->Fill(nCls);
+  
   delete clusterArray;
 
 }
@@ -779,7 +983,7 @@ void AliTRDQADataMakerRec::StartOfDetectorCycle()
 }
 
 //__________________________________________________________________________
-Int_t AliTRDQADataMakerRec::CheckPointer(TObject *obj, const char *name) 
+Int_t AliTRDQADataMakerRec::CheckPointer(TObject * const obj, const char *name) 
 {
   //
   // Checks initialization of pointers
@@ -789,7 +993,8 @@ Int_t AliTRDQADataMakerRec::CheckPointer(TObject *obj, const char *name)
   return !!obj;
 }
 //__________________________________________________________________________
-void AliTRDQADataMakerRec::BuildRatio(TH1D *ratio, TH1D *histN, TH1D*histD) {
+void AliTRDQADataMakerRec::BuildRatio(TH1D *ratio, TH1D * const histN, TH1D * const histD) 
+{
   //
   // Calculate the ratio of two histograms 
   // error are calculated assuming the histos have the same counts
@@ -820,4 +1025,25 @@ void AliTRDQADataMakerRec::BuildRatio(TH1D *ratio, TH1D *histN, TH1D*histD) {
   ratio->SetMaximum(1.1);
   ratio->SetMarkerStyle(20);
 }
+
+//__________________________________________________________________________
+Int_t AliTRDQADataMakerRec::FillBits(TH1D *hist, Int_t code, Int_t offset) 
+{
+  //
+  // Fill bits
+  //
+
+  Int_t nb = 0;
+  UInt_t test = 1;
+  for(Int_t i=0; i<8; i++) {
+    if (code & test) {
+      hist->Fill(i+offset);
+      nb++;
+    }
+    test *= 2;       
+  }
+  
+  return nb;
+}
+
 //__________________________________________________________________________