+
+//_____________________________________________________________________________
+Int_t AliTRDanalyzeHits()
+{
+ //
+ // Analyzes the hits and fills QA-histograms
+ //
+
+ Int_t rc = 0;
+
+ if (!gAlice) {
+ cout << "<AliTRDanalyzeHits> No AliRun object found" << endl;
+ rc = 1;
+ return rc;
+ }
+ gAlice->GetEvent(0);
+
+ // Get the pointer to the TRD detector
+ AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
+ if (!trd) {
+ cout << "<AliTRDanalyzeHits> No TRD detector found" << endl;
+ rc = 2;
+ return rc;
+ }
+
+ // Get the pointer to the geometry object
+ AliTRDgeometryFull *geo;
+ if (trd) {
+ geo = (AliTRDgeometryFull *) trd->GetGeometry();
+ }
+ else {
+ cout << "<AliTRDanalyzeHits> No TRD geometry found" << endl;
+ rc = 3;
+ return rc;
+ }
+
+ AliTRDparameter *par = new AliTRDparameter("TRDparameter","TRD parameter class");
+
+ // Define the histograms
+ TH1F *hQdedx = new TH1F("hQdedx","Charge dedx-hits",100,0.0,1000.0);
+ TH1F *hQtr = new TH1F("hQtr" ,"Charge TR-hits" ,100,0.0,1000.0);
+
+ Float_t rmin = geo->Rmin();
+ Float_t rmax = geo->Rmax();
+ Float_t length = geo->GetChamberLength(0,2);
+ Float_t width = geo->GetChamberWidth(0);
+ Int_t ncol = par->GetColMax(0);
+ Int_t nrow = par->GetRowMax(0,2,13);
+ Int_t ntime = ((Int_t) (rmax - rmin) / par->GetTimeBinSize());
+
+ TH2F *hZY = new TH2F("hZY" ,"Y vs Z (chamber 0)", nrow,-length/2.,length/2.
+ ,ntime, rmin, rmax);
+ TH2F *hXZ = new TH2F("hXZ" ,"Z vs X (plane 0)" , ncol, -width/2., width/2.
+ , nrow,-length/2.,length/2.);
+
+ TH1F *hNprimP = new TH1F("hNprimP","Number of primary electrons per cm (MIP pion)"
+ ,50,0.0,100.0);
+ TH1F *hNprimE = new TH1F("hNprimE","Number of primary electrons per cm (3GeV electron)"
+ ,50,0.0,100.0);
+ TH1F *hNtotP = new TH1F("hNtotP" ,"Number of electrons per cm (MIP pion)"
+ ,50,0.0,1000.0);
+ TH1F *hNtotE = new TH1F("hNtotE" ,"Number of electrons per cm (3GeV electron)"
+ ,50,0.0,1000.0);
+
+ // Get the pointer hit tree
+ TTree *hitTree = gAlice->TreeH();
+ if (!hitTree) {
+ cout << "<AliTRDanalyzeHits> No hit tree found" << endl;
+ rc = 4;
+ return rc;
+ }
+
+ Int_t countHits = 0;
+ Int_t nBytes = 0;
+
+ // Get the number of entries in the hit tree
+ // (Number of primary particles creating a hit somewhere)
+ Int_t nTrack = (Int_t) hitTree->GetEntries();
+ cout << "<AliTRDanalyzeHits> Found " << nTrack
+ << " primary particles with hits" << endl;
+
+ // Loop through all entries in the tree
+ for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
+
+ gAlice->ResetHits();
+ nBytes += hitTree->GetEvent(iTrack);
+
+ Int_t nPrimE = 0;
+ Int_t nPrimP = 0;
+ Int_t nTotE = 0;
+ Int_t nTotP = 0;
+
+ // Loop through the TRD hits
+ Int_t iHit = 0;
+ AliTRDhit *hit = (AliTRDhit *) trd->FirstHit(-1);
+ while (hit) {
+
+ countHits++;
+ iHit++;
+
+ Float_t x = hit->X();
+ Float_t y = hit->Y();
+ Float_t z = hit->Z();
+ Float_t q = hit->GetCharge();
+ Int_t track = hit->Track();
+ Int_t det = hit->GetDetector();
+ Int_t plane = geo->GetPlane(det);
+
+ if (q > 0) {
+ hQdedx->Fill(q);
+ hZY->Fill(z,y);
+ if (plane == 0) {
+ hXZ->Fill(x,z);
+ }
+ }
+ else if (q < 0) {
+ hQtr->Fill(TMath::Abs(q));
+ }
+
+ TParticle *part = gAlice->Particle(track);
+
+ if ((plane == 0) && (q > 0)) {
+
+ // 3 GeV electrons
+ if ((part->GetPdgCode() == 11) &&
+ (part->P() > 2.9)) {
+ nPrimE++;
+ nTotE += ((Int_t) q);
+ }
+
+ // MIP pions
+ if ((part->GetPdgCode() == -211) &&
+ (part->P() > 0.5)) {
+ nPrimP++;
+ nTotP += ((Int_t) q);
+ }
+
+ }
+
+ hit = (AliTRDhit *) trd->NextHit();
+
+ }
+
+ if (nPrimE > 0) hNprimE->Fill(((Double_t) nPrimE)/3.7);
+ if (nPrimP > 0) hNprimP->Fill(((Double_t) nPrimP)/3.7);
+ if (nTotE > 0) hNtotE->Fill(((Double_t) nTotE)/3.7);
+ if (nTotP > 0) hNtotP->Fill(((Double_t) nTotP)/3.7);
+
+ }
+
+ cout << "<AliTRDanalyzeHits> Found " << countHits << " hits in total" << endl;
+
+ TCanvas *cHits = new TCanvas("cHits","AliTRDanalyzeHits",50,50,600,600);
+ cHits->Divide(2,2);
+ cHits->cd(1);
+ hXZ->Draw("COL");
+ cHits->cd(2);
+ hZY->Draw("COL");
+ cHits->cd(3);
+ gPad->SetLogy();
+ hQdedx->Draw();
+ cHits->cd(4);
+ gPad->SetLogy();
+ hQtr->Draw();
+
+ TCanvas *cNel = new TCanvas("cNel","Ionization",50,50,600,600);
+ cNel->Divide(2,2);
+ cNel->cd(1);
+ hNprimE->SetStats();
+ hNprimE->Draw();
+ cNel->cd(2);
+ hNprimP->SetStats();
+ hNprimP->Draw();
+ cNel->cd(3);
+ hNtotE->SetStats();
+ hNtotE->Draw();
+ cNel->cd(4);
+ hNtotP->SetStats();
+ hNtotP->Draw();
+
+ TFile *fout = new TFile("TRD_ionization.root","RECREATE");
+ hNprimE->Write();
+ hNprimP->Write();
+ hNtotE->Write();
+ hNtotP->Write();
+ fout->Close();
+
+ return rc;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDcreateDigits()
+{
+ //
+ // Creates the digits from the hits of the slow simulator
+ //
+
+ Int_t rc = 0;
+
+ if (!gAlice) {
+ cout << "<AliTRDcreateDigits> No AliRun object found" << endl;
+ rc = 1;
+ return rc;
+ }
+ gAlice->GetEvent(0);
+
+ // Create the TRD digitzer
+ AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
+ ,"TRD digitizer class");
+ digitizer->InitDetector();
+
+ // Set the parameter
+ digitizer->SetDebug(1);
+
+ // Define the parameter object
+ // If no external parameter object is defined,
+ // default parameter will be used
+ AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
+ ,"TRD parameter class");
+ digitizer->SetParameter(parameter);
+
+ // Create the digits
+ if (!(digitizer->MakeDigits())) {
+ rc = 2;
+ return rc;
+ }
+
+ // Write the digits into the input file
+ if (!(digitizer->MakeBranch())) {
+ rc = 3;
+ return rc;
+ }
+
+ // Write the digits into the input file
+ if (!(digitizer->WriteDigits())) {
+ rc = 4;
+ return rc;
+ }
+
+ // Save the parameter object in the AliROOT file
+ if (!(parameter->Write())) {
+ rc = 4;
+ return rc;
+ }
+
+ return rc;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDanalyzeDigits()
+{
+ //
+ // Analyzes the digits
+ //
+
+ Int_t rc = 0;
+
+ const Int_t kNpla = AliTRDgeometry::Nplan();
+
+ if (!gAlice) {
+ cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
+ rc = 1;
+ return rc;
+ }
+ Int_t nPart = gAlice->GetEvent(0);
+
+ // Get the pointer to the TRD detector
+ AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
+ if (!trd) {
+ cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl;
+ rc = 2;
+ return rc;
+ }
+
+ // Get the digitizer object
+ TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
+ AliTRDparameter *parameter = (AliTRDparameter *) file->Get("TRDparameter");
+ if (!parameter) {
+ cout << "<AliTRDanalyzeDigits> No parameter object found" << endl;
+ rc = 3;
+ return rc;
+ }
+
+ // Define the histograms
+ Int_t adcRange = ((Int_t) parameter->GetADCoutRange());
+ TH1F *hAmpAll = new TH1F("hAmpAll" ,"Amplitude of the digits (all)"
+ ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
+ TH1F *hAmpEl = new TH1F("hAmpEl" ,"Amplitude of the digits (electrons)"
+ ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
+ TH1F *hAmpPi = new TH1F("hAmpPi" ,"Amplitude of the digits (pions)"
+ ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
+ TH1F *hAmpNoise = new TH1F("hAmpNoise","Amplitude of the digits (noise)"
+ ,5,-0.5,4.5);
+
+ // Get the pointer to the geometry object
+ AliTRDgeometry *geo;
+ if (trd) {
+ geo = trd->GetGeometry();
+ }
+ else {
+ cout << "<AliTRDanalyzeDigits> No TRD geometry found" << endl;
+ rc = 4;
+ return rc;
+ }
+
+ // Create the digits manager
+ AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
+
+ digitsManager->Open("TRD_test.root");
+
+ // Read the digits from the file
+ if (!(digitsManager->ReadDigits())) {
+ cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl;
+ rc = 5;
+ return rc;
+ }
+
+ AliTRDmatrix *matrix;
+
+ Int_t countDigits = 0;
+ Int_t iSec = trd->GetSensSector();
+ Int_t iCha = trd->GetSensChamber();
+ Int_t timeMax = geo->GetTimeTotal();
+
+ TProfile *hAmpTimeEl = new TProfile("hAmpTimeEl","Amplitude of the digits (electrons)"
+ ,timeMax,-0.5,((Double_t) timeMax)-0.5);
+ TProfile *hAmpTimePi = new TProfile("hAmpTimePi","Amplitude of the digits (pions)"
+ ,timeMax,-0.5,((Double_t) timeMax)-0.5);
+
+ // Loop over all planes
+ for (Int_t iPla = 0; iPla < kNpla; iPla++) {
+
+ Int_t iDet = geo->GetDetector(iPla,iCha,iSec);
+ Int_t rowMax = parameter->GetRowMax(iPla,iCha,iSec);
+ Int_t colMax = parameter->GetColMax(iPla);
+
+ if (iPla == 0) {
+ matrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
+ }
+
+ // Loop through the detector pixel
+ for (Int_t time = 0; time < timeMax; time++) {
+ for (Int_t col = 0; col < colMax; col++) {
+ for (Int_t row = 0; row < rowMax; row++) {
+
+ AliTRDdigit *digit = digitsManager->GetDigit(row,col,time,iDet);
+ Int_t amp = digit->GetAmp();
+ Int_t track0 = digitsManager->GetTrack(0,row,col,time,iDet);
+ Int_t track1 = digitsManager->GetTrack(1,row,col,time,iDet);
+ TParticle *particle = 0;
+ if (track0 > -1) {
+ particle = gAlice->Particle(track0);
+ }
+
+ if (amp > 0) {
+ countDigits++;
+ if (iPla == 0) {
+ matrix->SetSignal(row,col,time,amp);
+ }
+ }
+
+ // Total spectrum
+ hAmpAll->Fill(amp);
+
+ // Noise spectrum
+ if (track0 < 0) {
+ hAmpNoise->Fill(amp);
+ }
+
+ // Electron digit
+ if ((particle) && (particle->GetPdgCode() == 11) && (track1 < 0)) {
+ hAmpEl->Fill(amp);
+ hAmpTimeEl->Fill(time,amp);
+ }
+
+ // Pion digit
+ if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
+ hAmpPi->Fill(amp);
+ hAmpTimePi->Fill(time,amp);
+ }
+
+ delete digit;
+
+ }
+ }
+ }
+
+ }
+
+ cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl;
+
+ // Display the detector matrix
+ matrix->Draw();
+
+ TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",50,50,600,800);
+ cDigits->Divide(2,3);
+ cDigits->cd(1);
+ gPad->SetLogy();
+ hAmpAll->SetXTitle("Amplitude (ADC-channels)");
+ hAmpAll->SetYTitle("Entries");
+ hAmpAll->Draw();
+ cDigits->cd(2);
+ gPad->SetLogy();
+ hAmpNoise->SetXTitle("Amplitude (ADC-channels)");
+ hAmpNoise->SetYTitle("Entries");
+ hAmpNoise->Draw();
+ cDigits->cd(3);
+ gPad->SetLogy();
+ hAmpEl->SetXTitle("Amplitude (ADC-channels)");
+ hAmpEl->SetYTitle("Entries");
+ hAmpEl->Draw();
+ cDigits->cd(4);
+ gPad->SetLogy();
+ hAmpPi->SetXTitle("Amplitude (ADC-channels)");
+ hAmpPi->SetYTitle("Entries");
+ hAmpPi->Draw();
+ cDigits->cd(5);
+ hAmpTimeEl->SetXTitle("Timebin number");
+ hAmpTimeEl->SetYTitle("Mean amplitude");
+ hAmpTimeEl->Draw("HIST");
+ cDigits->cd(6);
+ hAmpTimePi->SetXTitle("Timebin number");
+ hAmpTimePi->SetYTitle("Mean amplitude");
+ hAmpTimePi->Draw("HIST");
+
+ TFile *fileOut = new TFile("digits_test.root","RECREATE");
+ hAmpAll->Write();
+ hAmpNoise->Write();
+ hAmpEl->Write();
+ hAmpPi->Write();
+ hAmpTimeEl->Write();
+ hAmpTimePi->Write();
+ fileOut->Close();
+
+ return rc;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDcreateCluster()
+{
+ //
+ // Creates the cluster from the digits
+ //
+
+ Int_t rc = 0;
+
+ // Create the clusterizer
+ AliTRDclusterizerV1 *clusterizer = new AliTRDclusterizerV1("TRDclusterizer"
+ ,"TRD clusterizer class");
+ clusterizer->SetVerbose(1);
+
+ // Define the parameter object
+ // If no external parameter object is defined,
+ // default parameter will be used
+ AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
+ ,"TRD parameter class");
+ parameter->SetClusMaxThresh(0);
+ parameter->SetClusSigThresh(0);
+ clusterizer->SetParameter(parameter);
+
+ // Open the file
+ if (!(clusterizer->Open("TRD_test.root",0))) {
+ rc = 1;
+ return rc;
+ }
+
+ // Load the digits
+ if (!(clusterizer->ReadDigits())) {
+ rc = 2;
+ return rc;
+ }
+
+ // Find the cluster
+ if (!(clusterizer->MakeClusters())) {
+ rc = 3;
+ return rc;
+ }
+
+ // Write the cluster tree into the file
+ if (!(clusterizer->WriteClusters(-1))) {
+ rc = 4;
+ return rc;
+ }
+
+ // Save the clusterizer class in the file
+ if (!(clusterizer->Write())) {
+ rc = 5;
+ return rc;
+ }
+
+ return rc;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDanalyzeCluster()
+{
+ //
+ // Analyzes the cluster
+ //
+
+ Int_t rc = 0;
+
+ if (!gAlice) {
+ cout << "<AliTRDanalyzeCluster> No AliRun object found" << endl;
+ rc = 1;
+ return rc;
+ }
+ gAlice->GetEvent(0);
+
+ // Get the pointer to the TRD detector
+ AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
+ if (!trd) {
+ cout << "<AliTRDanalyzeCluster> No TRD detector found" << endl;
+ rc = 2;
+ return rc;
+ }
+
+ // Define the histograms
+ TH1F *hClusAll = new TH1F("hClusAll" ,"Amplitude of the cluster (all)"
+ ,501,-0.5,500.5);
+ TH1F *hClusNoise = new TH1F("hClusNoise","Amplitude of the cluster (noise)"
+ , 5,-0.5, 4.5);
+ TH1F *hClusEl = new TH1F("hClusEl" ,"Amplitude of the cluster (electron)"
+ ,501,-0.5,500.5);
+ TH1F *hClusPi = new TH1F("hClusPi" ,"Amplitude of the cluster (pion)"
+ ,501,-0.5,500.5);
+
+ // Get the pointer to the geometry object
+ AliTRDgeometry *geo;
+ if (trd) {
+ geo = trd->GetGeometry();
+ }
+ else {
+ cout << "<AliTRDanalyzeCluster> No TRD geometry found" << endl;
+ rc = 3;
+ return rc;
+ }
+
+ // Get the pointer to the hit-tree
+ TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
+ TTree *clusterTree = (TTree *) file->Get("TreeR0_TRD");
+ if (!(clusterTree)) {
+ cout << "<AliTRDanalyzeCluster> No tree with clusters found" << endl;
+ rc = 4;
+ return rc;
+ }
+
+ // Get the pointer to the hit container
+ TObjArray *clusterArray = trd->RecPoints();
+ if (!(clusterArray)) {
+ cout << "<AliTRDanalyzeCluster> No clusterArray found" << endl;
+ rc = 5;
+ return rc;
+ }
+
+ // Set the branch address
+ clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
+ Int_t nEntries = clusterTree->GetEntries();
+ cout << "<AliTRDanalyzeCluster> Number of entries in the cluster tree = "
+ << nEntries
+ << endl;
+
+ Int_t countCluster = 0;
+ Int_t countOverlap = 0;
+
+ // Loop through all entries in the tree
+ Int_t nbytes;
+ for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
+
+ // Import the tree
+ nbytes += clusterTree->GetEvent(iEntry);
+
+ // Get the number of points in the detector
+ Int_t nCluster = clusterArray->GetEntriesFast();
+
+ // Loop through all TRD digits
+ for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
+
+ // Get the information for this digit
+ AliTRDcluster *cluster = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
+ Int_t detector = cluster->GetDetector();
+ Int_t sector = geo->GetSector(detector);
+ Int_t plane = geo->GetPlane(detector);
+ Int_t chamber = geo->GetChamber(detector);
+ Float_t energy = cluster->GetQ();
+ Int_t track0 = cluster->GetTrackIndex(0);
+ Int_t track1 = cluster->GetTrackIndex(1);
+ Int_t track2 = cluster->GetTrackIndex(2);
+ TParticle *particle = 0;
+ if (track0 > -1) {
+ particle = gAlice->Particle(track0);
+ }
+
+ countCluster++;
+ if (!cluster->Isolated()) countOverlap++;
+
+ // Total spectrum
+ hClusAll->Fill(energy);
+
+ if (cluster->Isolated()) {
+
+ // Noise spectrum
+ if (track0 < 0) {
+ hClusNoise->Fill(energy);
+ }
+
+ // Electron cluster
+ if ((particle) && (particle->GetPdgCode() == 11) && (track1 < 0)) {
+ hClusEl->Fill(energy);
+ }
+
+ // Pion cluster
+ if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
+ hClusPi->Fill(energy);
+ }
+
+ }
+
+ }
+
+ }
+
+ cout << "<AliTRDanalyzeCluster> Found " << countCluster << " cluster in total" << endl;
+ cout << "<AliTRDanalyzeCluster> Found " << countOverlap << " overlapping cluster" << endl;
+ cout << endl;
+
+ TCanvas *cCluster = new TCanvas("cCluster","AliTRDanalyzeCluster",50,50,600,600);
+ cCluster->Divide(2,2);
+
+ TF1 *fun;
+ cCluster->cd(1);
+ gPad->SetLogy();
+ hClusAll->Fit("landau","0");
+ fun = (TF1 *) hClusAll->GetListOfFunctions()->First();
+ Float_t meanAll = fun->GetParameter(1);
+ hClusAll->Draw();
+ fun->SetLineColor(2);
+ fun->Draw("SAME");
+
+ //cCluster->cd(2);
+ //gPad->SetLogy();
+ Float_t meanNoise = hClusNoise->GetMean();
+ //hClusNoise->Draw();
+
+ cCluster->cd(3);
+ gPad->SetLogy();
+ hClusEl->Fit("landau","0");
+ fun = (TF1 *) hClusEl->GetListOfFunctions()->First();
+ fun->SetLineColor(2);
+ Float_t meanEl = fun->GetParameter(1);
+ hClusEl->Draw();
+ fun->Draw("SAME");
+
+ cCluster->cd(4);
+ gPad->SetLogy();
+ hClusPi->Fit("landau","0");
+ fun = (TF1 *) hClusPi->GetListOfFunctions()->First();
+ fun->SetLineColor(2);
+ Float_t meanPi = fun->GetParameter(1);
+ hClusPi->Draw();
+ fun->Draw("SAME");
+
+ cout << endl;
+ cout << "##################################################################" << endl;
+ cout << " Mean all = " << meanAll << endl;
+ cout << " Mean noise = " << meanNoise << endl;
+ cout << " Mean electrons = " << meanEl << endl;
+ cout << " Mean pions = " << meanPi << endl;
+ cout << "##################################################################" << endl;
+ cout << endl;
+
+ return rc;
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff