X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDQADataMakerRec.cxx;h=89c5e57c6329e41a3eb3f4f399cb8e47c6b82276;hb=5c5d503aeb880dc107412dd3105750a3ac29ea26;hp=2b95b0565bc46997e1dc780c0d13b46c12d2770d;hpb=9b99c029c45a22703634c1d502945ec073e265e7;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDQADataMakerRec.cxx b/TRD/AliTRDQADataMakerRec.cxx index 2b95b0565bc..89c5e57c632 100644 --- a/TRD/AliTRDQADataMakerRec.cxx +++ b/TRD/AliTRDQADataMakerRec.cxx @@ -28,9 +28,9 @@ // --- ROOT system --- #include #include -#include -#include -#include +#include +#include +#include #include #include #include @@ -38,11 +38,15 @@ // --- AliRoot header files --- #include "AliESDEvent.h" #include "AliLog.h" +#include "AliRawReader.h" #include "AliTRDcluster.h" #include "AliTRDQADataMakerRec.h" #include "AliTRDgeometry.h" -#include "AliTRDdataArrayI.h" -#include "AliTRDrawStreamTB.h" +#include "AliTRDrawStream.h" + +#include "AliTRDdigitsManager.h" +#include "AliTRDSignalIndex.h" +#include "AliTRDarrayADC.h" #include "AliQAChecker.h" @@ -50,7 +54,7 @@ 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 +87,187 @@ 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(); - - AliInfo("End of TRD cycle"); - - if (task == AliQA::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); - } - - - // Rec points full chambers - for (Int_t i=0; i<540; i++) { - - //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; bGetXaxis()->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); + ResetEventTrigClasses(); + AliDebug(AliQAv1::GetQADebugLevel(), "End of TRD cycle"); + // + for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) { // loop over event types + // + if (!IsValidEventSpecie(specie, list)) continue; + SetEventSpecie(AliRecoParam::ConvertIndex(specie)); + // + for (int itc=-1;itcDivide(mnCls); } - - //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; bGetXaxis()->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; + // + if (task == AliQAv1::kRECPOINTS) { + // + TH1F * hist = new TH1F("fitHist", "", 200, -0.5, 199.5); + TObjArray& arrRP = *GetRecPointsDataOfTrigClass(itc); // RS Histos matching to trigger class + // fill detector map; + TH1* h0 = (TH1*) arrRP[0]; + TH1* detMap = (TH1*) arrRP[87]; + if (h0 && detMap) { + for(Int_t i = 0 ; i < 540 ; i++) { + Double_t v = h0->GetBinContent(i+1); + Int_t sm = i/30; + Int_t det = i%30; + Int_t bin = detMap->FindBin(sm, det); + detMap->SetBinContent(bin, v); } - //printf("v = %f\n", value); - hist->SetBinContent(b, svalue); } - - 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"); - - 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))); - } - } - - - // 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)); + // + // Rec points full chambers + TH2* h2tmp = (TH2*) arrRP[1]; + TH1* h12 = (TH1*) arrRP[12]; + TH1* h13 = (TH1*) arrRP[13]; + if (h2tmp && h12 && h13) { + for (Int_t i = 0 ; i < 540 ; i++) { + //AliDebug(AliQAv1::GetQADebugLevel(), Form("I = %d", i)); + hist->Reset(); + // project TH2F into TH1F + for(Int_t b = 1 ; b < hist->GetXaxis()->GetNbins()-1 ; b++) { + Double_t xvalue = hist->GetBinCenter(b); + Int_t bin = h2tmp->FindBin(i,xvalue); + Double_t value = h2tmp->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"); + h12->Fill(fit->GetParameter(1)); + h13->Fill(fit->GetParameter(2)); + } + } + // + // time-bin by time-bin sm by sm + TH3* h3tmp = (TH3*) arrRP[10]; + if (h3tmp) { + for(Int_t i=0; i<18; i++) { // loop over super-modules + for(Int_t j=0; jReset(); + 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 = h3tmp->FindBin(det,j,xvalue); + Double_t value = h3tmp->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"); + // + TH1 *hi = (TH1*)arrRP[14+18+i]; + if (!hi) continue; + hi->SetMarkerStyle(20); + Int_t bin = hi->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 + // + hi->SetBinContent(bin, value); + hi->SetBinError(bin, error); + } // j + } // i + // + // for numerical convergence + TF1 *form = new TF1("formLandau", "landau", 0, 200); + // + // time-bin by time-bin chamber by chamber + for (Int_t i=0; i<540; i++) { + for(Int_t j=0; jReset(); + for(Int_t b = 1 ; b < hist->GetXaxis()->GetNbins()-1 ; b++) { + Double_t xvalue = hist->GetBinCenter(b); + Int_t bin = h3tmp->FindBin(i,j,xvalue); + Double_t value = h3tmp->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; + TH1 *hi = (TH1*)arrRP[14+sm]; + if (!hi) continue; + Int_t bin = hi->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; + // + hi->SetBinContent(bin, value); + hi->SetBinError(bin, error); + } // j + } // i + } // h3tmp + if (hist) delete hist; + } // RECPOINTS + // + ////////////////////////// + // const Int_t knbits = 6; + // const char *suf[knbits] = {"TPCi", "TPCo", "TPCz", "TRDo", "TRDr", "TRDz"}; + //const char *sufRatio[4] = {"TRDrTRDo", "TRDoTPCo", "TRDrTPCo", "TRDzTPCo"}; - 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; bGetXaxis()->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); + 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}; + // + TObjArray& arrES = *GetESDsDataOfTrigClass(itc); // RS Histos matching to trigger class + // create ratios + for(Int_t type = 0 ; type < 2 ; type++) { + for(Int_t i = 0 ; i < knRatio ; i++) { + TH1 *ratio = (TH1*)arrES[19 + 2*i + type]; + TH1 *histN = (TH1*)arrES[3 + 2*kN[i] + type]; + TH1 *histD = (TH1*)arrES[3 + 2*kD[i] + type]; + if ( ! (ratio && histN && histD) ) continue; + BuildRatio(ratio, histN, histD); + //ratio->Reset(); + //ratio->Add(histN); + //ratio->Divide(histD); + } } - - 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"); - - 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)); - } - } - - 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) { - - 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; iAt(19 + 2*i + type); - TH1D *histN = (TH1D*)list->At(3 + 2*kN[i] + type); - TH1D *histD = (TH1D*)list->At(3 + 2*kD[i] + type); - + // ratio for the fraction of electrons per stack + TH1 *histN = (TH1*)arrES[33]; + TH1 *histD = (TH1*)arrES[32]; + TH1 *ratio = (TH1*)arrES[34]; BuildRatio(ratio, histN, histD); - //ratio->Reset(); - //ratio->Add(histN); - //ratio->Divide(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); - } - - + } // ESDs + } // RS: loop over eventual clones per trigger class + } // loop over species // call the checker - AliQAChecker::Instance()->Run(AliQA::kTRD, task, list) ; + AliQAChecker::Instance()->Run(AliQAv1::kTRD, task, list) ; + // } //____________________________________________________________________________ @@ -250,70 +276,96 @@ void AliTRDQADataMakerRec::InitESDs() // // Create ESDs histograms in ESDs subdir // - - const Int_t kNhist = 36+5; + 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 TH1F("qaTRD_esd_ntracks", "TRD esd ntracks;Number of tracks;Counts", 300, -0.5, 299.5); + hist[++histoCounter] = new TH1F("qaTRD_esd_sector", "TRD esd sector;Sector;Counts", 18, -0.5, 17.7); + hist[++histoCounter] = new TH1F("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; iSumw2(); - Add2ESDsList(hist[i], i); + Add2ESDsList(hist[i], i, !expert, image); } - + // + ClonePerTrigClass(AliQAv1::kESDS); // this should be the last line } //____________________________________________________________________________ @@ -322,70 +374,87 @@ 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 TH1F("qaTRD_recPoints_det", "RRD recPoints det;Detector ID of the cluster;Counts", 540, -0.5, 539.5); + hist[1] = new TH2F("qaTRD_recPoints_amp", "TRD recPoints amp;Amplitude;??", 540, -0.5, 539, 200, -0.5, 199.5); + hist[2] = new TH1F("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 TH1F("qaTRD_recPoints_dist2", "TRD recPoints dist2;residuals [2pad];Counts", 100, -1, 1); + hist[4] = new TH1F("qaTRD_recPoints_dist3", "TRD recPoints dist3;residuals [3pad];Counts", 100, -1, 1); + hist[5] = new TH1F("qaTRD_recPoints_dist4", "TRD recPoints dist4;residuals [4pad];Counts", 100, -1, 1); + hist[6] = new TH1F("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 TH2F("qaTRD_recPoints_rowCol", "TRDrecPointsrowCol;row;col", 16, -0.5, 15.5, 145, -0.5, 144.5); + hist[8] = new TH1F("qaTRD_recPoints_time", "TRDrecPoints time;time bin;Counts", kTimeBin, -0.5, kTimeBin-0.5); + hist[9] = new TH1F("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 TH3F("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 TH1F("qaTRD_recPoints_ampMPV", "TRD recPoints ampMPV;amplitude MPV;Counts", 150, 0, 150); + hist[13] = new TH1F("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); + hist[14+i] = new TH2F(Form("qaTRD_recPoints_sigTime_sm%d",i), Form("sm%d;det;time bin",i), + 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); + hist[14+18+i] = new TH1F(Form("qaTRD_recPoints_sigTimeShape_sm%d", i), + Form("sm%d;time bin;signal",i), + 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), + hist[50+i] = new TH1F(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), + hist[68+i] = new TH1F(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 TH1F("qaTRD_recPoints_signal", "TRD recPoints signal;amplitude;Counts", 400, -0.5, 399.5); + hist[87] = new TH2F("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 TH1F(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 TH2F(Form("qaTRD_recPoints_map%d", i), ";col;row", 16, -0.5, 15.5, 144, -0.5, 143.5); + //} for(Int_t i=0; iSumw2(); - Add2RecPointsList(hist[i], i); + Add2RecPointsList(hist[i], i, !expert, image); } + // + ClonePerTrigClass(AliQAv1::kRECPOINTS); // this should be the last line } //____________________________________________________________________________ @@ -394,30 +463,46 @@ 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 kNhist = 7; + TH1 *hist[kNhist]; - const Int_t kSM = 18; - //const Int_t kNCh = 540; - const Int_t kNhist = 4+kSM; - TH1D *hist[kNhist]; - - // 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); - hist[2] = new TH1D("qaTRD_raws_timeBin", ";time bin", 40, -0.5, 39.5); - hist[3] = new TH1D("qaTRD_raws_smId", ";supermodule", 18, -0.5, 17.5); - // + hist[0] = new TH2F("qaTRD_raws_nADC","number of ADC channels;sector;detector", 18, -0.5, 17.5, 30, -0.5, 29.5); + hist[1] = new TH2F("qaTRD_raws_nCls", "number of clusters;sector;detector", 18, -0.5, 17.5, 30, -0.5, 29.5); + hist[2] = new TH2F("qaTRD_raws_meanSig", "mean signal;sector;detector", 18, -0.5, 17.5, 30, -0.5, 29.5); - // one double per MCM (not published) - const Int_t kNMCM = 30 * 8 * 16; - for(Int_t i=0; iSumw2(); - Add2RawsList(hist[i], i); + Add2RawsList(hist[i], i, !expert, image, !saveCorr); } - + // + ClonePerTrigClass(AliQAv1::kRAWS); // this should be the last line } //____________________________________________________________________________ @@ -426,14 +511,14 @@ void AliTRDQADataMakerRec::MakeESDs(AliESDEvent * esd) // // Make QA data from ESDs // - + Int_t nTracks = esd->GetNumberOfTracks(); - GetESDsData(0)->Fill(nTracks); + FillESDsData(0,nTracks); // track loop - for (Int_t i=0; iGetTrack(i); + AliESDtrack *track = esd->GetTrack(iTrack); const AliExternalTrackParam *paramOut = track->GetOuterParam(); const AliExternalTrackParam *paramIn = track->GetInnerParam(); @@ -454,7 +539,7 @@ void AliTRDQADataMakerRec::MakeESDs(AliESDEvent * esd) UInt_t u = 1; UInt_t status = track->GetStatus(); for(Int_t bit=0; bit<32; bit++) - if (u<Fill(bit); + if (u<Fill(pt); - GetESDsData(2*b+4)->Fill(extZ); + FillESDsData(2*b+3,pt); + FillESDsData(2*b+4,extZ); } } // clusters for(Int_t b=0; b<3; b++) - if (bit[3+b]) GetESDsData(b+15)->Fill(track->GetTRDncls0()); + if (bit[3+b]) FillESDsData(b+15,track->GetTRDncls0()); // refitted only if (!bit[4]) continue; @@ -487,38 +572,48 @@ void AliTRDQADataMakerRec::MakeESDs(AliESDEvent * esd) //fBudget->Fill(track->GetTRDBudget()); //fSignal->Fill(track->GetTRDsignal()); - GetESDsData(1)->Fill(sector); - GetESDsData(18)->Fill(track->GetP(), track->GetTRDsignal()); + FillESDsData(1,sector); + FillESDsData(18,track->GetP(), track->GetTRDsignal()); - GetESDsData(27)->Fill(track->GetTRDQuality()); - GetESDsData(28)->Fill(track->GetTRDBudget()); - GetESDsData(29)->Fill(track->GetTRDchi2()); - GetESDsData(30)->Fill(track->GetTRDTimBin(0)); - GetESDsData(31)->Fill(track->GetTRDpidQuality()); + FillESDsData(27,track->GetTRDQuality()); + FillESDsData(28,track->GetTRDBudget()); + FillESDsData(29,track->GetTRDchi2()); + FillESDsData(30,track->GetTRDTimBin(0)); + FillESDsData(31,track->GetTRDntrackletsPID()); + + // dedx + for(Int_t k=0; k<4; ++k) { + Double_t dedx = 0; + for(Int_t j=0; j<6; j++) { + dedx += track->GetTRDslice(j, k-1); + } + FillESDsData(41+k,paramOut->GetP(), dedx/6.); + } + // probabilities if (status & AliESDtrack::kTRDpid) { - for(Int_t i=0; iFill(track->GetTRDpid(i)); + for(Int_t k=0; kGetTRDpid(k)); } // probabilities uniformity - if (track->GetTRDpidQuality() < 6) continue; - GetESDsData(35)->Fill(paramOut->GetZ()/paramOut->GetX()); + if (track->GetTRDntrackletsPID() < 6) continue; + FillESDsData(35,paramOut->GetZ()/paramOut->GetX()); Int_t idx = 5 * sector + stack; - GetESDsData(32)->Fill(idx); // all tracks + FillESDsData(32,idx); // all tracks if (track->GetTRDpid(AliPID::kElectron) > 0.9) - GetESDsData(33)->Fill(idx); // electrons only + FillESDsData(33,idx); // electrons only /* - 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", 300, 0, 100); - hist[30] = new TH1D("qaTRD_esd_timeBin", 7, -0.5, 6.5); - hist[31] = new TH1D("qaTRD_esd_pidQuality", 7, -0.5, 6.5); + hist[27] = new TH1F("qaTRD_esd_quality", ";quality", 120, 0, 12); + hist[28] = new TH1F("qaTRD_esd_budget", ";NN", 110, -1000, 100); + hist[29] = new TH1F("qaTRD_esd_chi2", ";chi2", 300, 0, 100); + hist[30] = new TH1F("qaTRD_esd_timeBin", 7, -0.5, 6.5); + hist[31] = new TH1F("qaTRD_esd_pidQuality", 7, -0.5, 6.5); */ /* @@ -556,7 +651,10 @@ void AliTRDQADataMakerRec::MakeESDs(AliESDEvent * esd) */ } - + // + IncEvCountCycleESDs(); + IncEvCountTotalESDs(); + // } //______________________________________________________________________________ @@ -615,47 +713,108 @@ void AliTRDQADataMakerRec::MakeRaws(AliRawReader* rawReader) // Makes QA data from raw data // - // 157 - // T9 -- T10 - - //const Int_t kSM = 18; - //const Int_t kROC = 30; - const Int_t kROB = 8; - //const Int_t kLayer = 6; - //const Int_t kStack = 5; - const Int_t kMCM = 16; - // const Int_t kADC = 22; + AliInfo("Making QA for Raws"); + + // take histograms (RS: arrays of histos) + static TObjArray mnADC,mnCls,mClsDet,mADC,mCls,mClsTb,mClsDetAmp; + GetMatchingRawsData(0,&mnADC); + GetMatchingRawsData(1,&mnCls); + GetMatchingRawsData(2,&mClsDet); + GetMatchingRawsData(3,&mADC); + GetMatchingRawsData(4,&mCls); + GetMatchingRawsData(5,&mClsTb); + GetMatchingRawsData(6,&mClsDetAmp); + + const Int_t baseline = 10; + + // configure the reader + rawReader->Reset(); + rawReader->SelectEquipment(0, 1024, 1041); + rawReader->Select("TRD"); + + AliTRDrawStream *data = new AliTRDrawStream(rawReader); + + // build data manager + AliTRDdigitsManager *digitsManager; + digitsManager = new AliTRDdigitsManager(kTRUE); + digitsManager->CreateArrays(); + + // error container + const UShort_t kErrorChmb = 1411; + UShort_t **fErrorContainer = new UShort_t *[2]; + fErrorContainer[0] = new UShort_t[kErrorChmb]; + fErrorContainer[1] = new UShort_t[kErrorChmb]; + + Int_t det = 0; + Int_t row, col; - AliTRDrawStreamTB *raw = new AliTRDrawStreamTB(rawReader); + while ((det = data->NextChamber(digitsManager, NULL, fErrorContainer)) >= 0){ + + //printf("DET = %d\n", det); - raw->SetRawVersion(3); - raw->Init(); + AliTRDSignalIndex* indexes = digitsManager->GetIndexes(det); + if (indexes->HasEntry()) { + + AliTRDarrayADC *digits = (AliTRDarrayADC*) digitsManager->GetDigits(det); - while (raw->Next()) { + while(indexes->NextRCIndex(row, col)) { - GetRawsData(0)->Fill(raw->GetDet()); + for (int ih=mnADC.GetEntriesFast();ih--;) ((TH2*)mnADC.UncheckedAt(ih))->Fill(det/30, det%30); + + for(Int_t tb = 0; tb < digits->GetNtime(); tb++) { + Int_t value = digits->GetData(row, col, tb); + for (int ih=mADC.GetEntriesFast();ih--;) ((TH1*)mADC.UncheckedAt(ih))->Fill(value); + + // simple clusterizer + if (col < 1 || col > digits->GetNcol()-2) continue; + if (tb < 1 || tb > digits->GetNtime()-2) continue; + + value -= baseline; - // possibly needs changes with the new reader !! - Int_t *sig = raw->GetSignals(); - for(Int_t i=0; i<3; i++) GetRawsData(1)->Fill(sig[i]); - // --- + Int_t valueL = digits->GetData(row, col-1, tb) - baseline; + if (valueL >= value) continue; + + Int_t valueR = digits->GetData(row, col+1, tb) - baseline; + if (valueR >= value) continue; + + Int_t valueUp = digits->GetData(row, col-1, tb+1) + + digits->GetData(row, col, tb+1) + digits->GetData(row, col+1, tb+1) - 3 * baseline; + if (valueUp < 10) continue; + + Int_t valueDown = digits->GetData(row, col-1, tb-1) + + digits->GetData(row, col, tb-1) + digits->GetData(row, col+1, tb-1) - 3 * baseline; + if (valueDown < 10) continue; + + Int_t valueTot = value + valueL + valueR; + if (valueTot < 0) continue; - GetRawsData(2)->Fill(raw->GetTimeBin()); + for (int ih=mCls.GetEntriesFast();ih--;) ((TH1*)mCls.UncheckedAt(ih))->Fill(valueTot); + for (int ih=mClsTb.GetEntriesFast();ih--;) ((TH2*)mClsTb.UncheckedAt(ih))->Fill(tb, valueTot); + for (int ih=mClsDetAmp.GetEntriesFast();ih--;) ((TH2*)mClsDetAmp.UncheckedAt(ih))->Fill(det, valueTot); - // calculate the index; - Int_t sm = raw->GetSM(); - Int_t roc = raw->GetROC(); - Int_t rob = raw->GetROB(); - Int_t mcm = raw->GetMCM(); - //Int_t adc = raw->GetADC(); + if (valueTot < 200) { + for (int ih=mnCls.GetEntriesFast();ih--;) ((TH2*)mnCls.UncheckedAt(ih))->Fill(det/30, det%30); + for (int ih=mClsDet.GetEntriesFast();ih--;) ((TH2*)mClsDet.UncheckedAt(ih))->Fill(det/30, det%30, valueTot); + } + + } + } - //Int_t index = roc * (kROB*kMCM*kADC) + rob * (kMCM*kADC) + mcm * kADC + adc; - Int_t index = roc * (kROB*kMCM) + rob * kMCM + mcm; - GetRawsData(3)->Fill(sm); - GetRawsData(4+sm)->Fill(index); + digitsManager->ClearArrays(det); // do we need this if object will be deleted ?? + } } - - delete raw; + + delete [] fErrorContainer[0]; + delete [] fErrorContainer[1]; + delete [] fErrorContainer; + fErrorContainer = NULL; + + delete digitsManager; + delete data; + // + IncEvCountCycleRaws(); + IncEvCountTotalRaws(); + // } //____________________________________________________________________________ @@ -666,7 +825,7 @@ void AliTRDQADataMakerRec::MakeRecPoints(TTree * clustersTree) // // Info("MakeRecPoints", "making"); - + Int_t nsize = Int_t(clustersTree->GetTotBytes() / (sizeof(AliTRDcluster))); TObjArray *clusterArray = new TObjArray(nsize+1000); @@ -678,42 +837,88 @@ 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; + } + */ + + ///// + TObjArray *hists3D = GetMatchingRecPointsData(10); //RS no alias for 3d histo filling, to directly + // 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(1+c->GetNPads())->Fill(c->GetCenter()); + FillRecPointsData(0,iDet); + FillRecPointsData(86,c->GetQ()); + FillRecPointsData(1,iDet, c->GetQ()); + FillRecPointsData(2,nPads); + if (nPads < 6) + FillRecPointsData(1+c->GetNPads(),c->GetCenter()); + + if (nPads < 10) + FillRecPointsData(88+nPads-2,c->GetQ()); + else FillRecPointsData(96,c->GetQ()); //if (c->GetPadTime() < 5) - ((TH2D*)GetRecPointsData(7))->Fill(c->GetPadRow(), c->GetPadCol()); - GetRecPointsData(8)->Fill(c->GetPadTime()); - - ((TH3D*)GetRecPointsData(10))->Fill(iDet, c->GetPadTime(), c->GetQ()); + FillRecPointsData(7,c->GetPadRow(), c->GetPadCol()); + FillRecPointsData(8,c->GetPadTime()); + // + if (hists3D) { + for (int ih=hists3D->GetEntriesFast();ih--;) { + TH3F * ahist = dynamic_cast(hists3D->At(ih)); + if (ahist) ahist->Fill(iDet, c->GetPadTime(), c->GetQ()); + } + } Int_t iSM = iDet / 30; - GetRecPointsData(50+iSM)->Fill(c->GetPadTime()); - GetRecPointsData(68+iSM)->Fill(c->GetPadTime(), c->GetQ()); + FillRecPointsData(50+iSM,c->GetPadTime()); + FillRecPointsData(68+iSM,c->GetPadTime(), c->GetQ()); + + // total charge sm / det / timeBin + // FillRecPointsData(14+iSM,iDet-iSM*30, c->GetPadTime(), c->GetQ()); + // PRF for 2pad //if (c->GetNPads() == 2) { @@ -726,17 +931,24 @@ void AliTRDQADataMakerRec::MakeRecPoints(TTree * clustersTree) if (sig[0] == 0 && sig[1] == 0 && sig[4] == 0 && sig[5] == 0 && sig[6] == 0) frac = -1. * sig[2] / (sig[2] + sig[3]); - if (frac > -10) ((TProfile*)GetRecPointsData(11))->Fill(c->GetCenter(), frac); + if (frac > -10) FillRecPointsData(11,c->GetCenter(), frac); //} } + clusterArray->Delete(); } + /* for(Int_t i=0; i<540; i++) - if (nDet[i] > 0) GetRecPointsData(9)->Fill(nDet[i]); - + if (nDet[i] > 0) FillRecPointsData(9,nDet[i]); + */ + FillRecPointsData(9,nCls); + delete clusterArray; - + // + IncEvCountCycleRecPoints(); + IncEvCountTotalRecPoints(); + // } //____________________________________________________________________________ @@ -759,7 +971,7 @@ Int_t AliTRDQADataMakerRec::CheckPointer(TObject *obj, const char *name) return !!obj; } //__________________________________________________________________________ -void AliTRDQADataMakerRec::BuildRatio(TH1D *ratio, TH1D *histN, TH1D*histD) { +void AliTRDQADataMakerRec::BuildRatio(TH1 *ratio, TH1 *histN, TH1* histD) { // // Calculate the ratio of two histograms // error are calculated assuming the histos have the same counts @@ -791,3 +1003,20 @@ void AliTRDQADataMakerRec::BuildRatio(TH1D *ratio, TH1D *histN, TH1D*histD) { ratio->SetMarkerStyle(20); } //__________________________________________________________________________ + +Int_t AliTRDQADataMakerRec::FillBits(TH1F *hist, Int_t code, Int_t offset) { + + 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; +} + +//__________________________________________________________________________