From fa148e6c745790262f4326f73145ce5d231f5fc1 Mon Sep 17 00:00:00 2001 From: cblume Date: Mon, 25 Mar 2002 20:03:21 +0000 Subject: [PATCH] Update of macros --- TRD/AliTRDanaCluster.C | 95 +++++ TRD/AliTRDanaDigits.C | 115 ++++++ TRD/AliTRDdigits2cluster.C | 47 +++ TRD/AliTRDdisplayDigits3D.C | 2 +- TRD/AliTRDdisplayTracks.C | 143 ++++++++ TRD/AliTRDhits2digits.C | 47 +++ TRD/AliTRDhits2sdigits.C | 49 +++ TRD/AliTRDmerge.C | 27 ++ TRD/AliTRDrunSimple.C | 274 +++++++++++++++ TRD/AliTRDsdigits2digits.C | 54 +++ TRD/AliTRDtest.C | 683 +++++++++++++++++++++++++++++++++++- 11 files changed, 1530 insertions(+), 6 deletions(-) create mode 100644 TRD/AliTRDanaCluster.C create mode 100644 TRD/AliTRDanaDigits.C create mode 100644 TRD/AliTRDdigits2cluster.C create mode 100644 TRD/AliTRDdisplayTracks.C create mode 100644 TRD/AliTRDhits2digits.C create mode 100644 TRD/AliTRDhits2sdigits.C create mode 100644 TRD/AliTRDmerge.C create mode 100644 TRD/AliTRDrunSimple.C create mode 100644 TRD/AliTRDsdigits2digits.C diff --git a/TRD/AliTRDanaCluster.C b/TRD/AliTRDanaCluster.C new file mode 100644 index 00000000000..a1b3315bfd0 --- /dev/null +++ b/TRD/AliTRDanaCluster.C @@ -0,0 +1,95 @@ +void AliTRDanaCluster() +{ + +///////////////////////////////////////////////////////////////////////// +// +// Example macro for the analysis of the TRD cluster +// +///////////////////////////////////////////////////////////////////////// + + // Dynamically link some shared libs + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + } + + // Input file name + Char_t *alifile = "galice.root"; + + // Event number + Int_t nEvent = 0; + + // Define the histograms + TH1F *hCharge = new TH1F("hCharge","Cluster charge",100,0.0,1000.0); + + // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits + TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile); + if (!gafl) { + cout << "Open the ALIROOT-file " << alifile << endl; + gafl = new TFile(alifile); + } + else { + cout << alifile << " is already open" << endl; + } + + // Get AliRun object from file or create it if not on file + gAlice = (AliRun*) gafl->Get("gAlice"); + if (gAlice) { + cout << "AliRun object found on file" << endl; + } + else { + gAlice = new AliRun("gAlice","Alice test program"); + } + + // Import the Trees for the event nEvent in the file + Int_t nparticles = gAlice->GetEvent(nEvent); + if (nparticles <= 0) break; + + // Get the pointer to the hit-tree + Char_t treeName[14]; + sprintf(treeName,"TreeR%d_TRD",nEvent); + TTree *clusterTree = gafl->Get(treeName); + clusterTree->Print(); + // Get the pointer to the detector classes + AliTRDv1 *trd = (AliTRDv1*) gAlice->GetDetector("TRD"); + // Get the geometry + AliTRDgeometry *geo = trd->GetGeometry(); + // Get the pointer to the hit container + TObjArray *clusterArray = trd->RecPoints(); + // Set the branch address + clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray); + + Int_t nEntries = clusterTree->GetEntries(); + cout << "nEntries = " << nEntries << endl; + + // 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 charge = cluster->GetQ(); + + hCharge->Fill(charge); + + } + + } + + TCanvas *c1 = new TCanvas("c1","Cluster",50,50,600,400); + hCharge->Draw(); + +} diff --git a/TRD/AliTRDanaDigits.C b/TRD/AliTRDanaDigits.C new file mode 100644 index 00000000000..7e3e9d7c428 --- /dev/null +++ b/TRD/AliTRDanaDigits.C @@ -0,0 +1,115 @@ +void AliTRDanaDigits() +{ + +///////////////////////////////////////////////////////////////////////// +// +// Example macro for the analysis of the TRD digits and the use +// of the AliTRDmatrix class. +// +///////////////////////////////////////////////////////////////////////// + + // Dynamically link some shared libs + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + } + + // Input file name + //Char_t *alifile = "galice.root"; + Char_t *alifile = "galice_jiri.root"; + + // Event number + Int_t nEvent = 0; + + // Define the objects + AliTRDv1 *trd; + AliTRDgeometry *geo; + AliTRDdigit *digit; + + Int_t track; + + // Connect the AliRoot file containing Geometry, Kine, Hits, and digits + TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile); + if (!gafl) { + cout << "Open the ALIROOT-file " << alifile << endl; + gafl = new TFile(alifile); + } + else { + cout << alifile << " is already open" << endl; + } + + // Get AliRun object from file or create it if not on file + gAlice = (AliRun*) gafl->Get("gAlice"); + if (gAlice) + cout << "AliRun object found on file" << endl; + else + gAlice = new AliRun("gAlice","Alice test program"); + + // Import the Trees for the event nEvent in the file + Int_t nparticles = gAlice->GetEvent(nEvent); + if (nparticles <= 0) break; + + // Get the pointer to the detector object + trd = (AliTRDv1*) gAlice->GetDetector("TRD"); + + // Get the pointer to the geometry object + if (trd) { + geo = trd->GetGeometry(); + } + else { + cout << "Cannot find the geometry" << endl; + break; + } + + // Create the digits manager + AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager(); + digitsManager->SetDebug(1); + digitsManager->SetSDigits(kTRUE); + + // Read the digits from the file + digitsManager->Open(alifile); + digitsManager->ReadDigits(); + + // Get the detector number + Int_t iDet = 423; + cout << " iDet = " << iDet << endl; + + // Define the detector matrix for one chamber + const Int_t iSec = geo->GetSector(iDet); + const Int_t iCha = geo->GetChamber(iDet); + const Int_t iPla = geo->GetPlane(iDet); + Int_t rowMax = geo->GetRowMax(iPla,iCha,iSec); + Int_t colMax = geo->GetColMax(iPla); + Int_t timeMax = geo->GetTimeMax(); + cout << "Geometry: rowMax = " << rowMax + << " colMax = " << colMax + << " timeMax = " << timeMax << endl; + AliTRDmatrix *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++) { + + digit = digitsManager->GetDigit(row,col,time,iDet); + track = digitsManager->GetTrack(0,row,col,time,iDet); + + printf("time=%d, col=%d, row=%d, amp=%f\n",time,col,row,digit->GetAmp()); + matrix->SetSignal(row,col,time,digit->GetAmp()); + + delete digit; + + } + } + } + + // Display the detector matrix + matrix->Draw(); + //matrix->DrawRow(18); + //matrix->DrawCol(58); + //matrix->DrawTime(20); + matrix->ProjRow(); + matrix->ProjCol(); + matrix->ProjTime(); + +} diff --git a/TRD/AliTRDdigits2cluster.C b/TRD/AliTRDdigits2cluster.C new file mode 100644 index 00000000000..cac782ea092 --- /dev/null +++ b/TRD/AliTRDdigits2cluster.C @@ -0,0 +1,47 @@ +void AliTRDdigits2cluster() +{ + +///////////////////////////////////////////////////////////////////////// +// +// Creates cluster from the digit information. +// +///////////////////////////////////////////////////////////////////////// + + // Dynamically link some shared libs + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + cout << "Loaded shared libraries" << endl; + } + + // Input and output file names + Char_t *infile = "galice.root"; + Char_t *outfile = "TRDclusters.root"; + + // Create the clusterizer + AliTRDclusterizerV1 *clusterizer = + new AliTRDclusterizerV1("clusterizer","Clusterizer class"); + + // Set the parameter + clusterizer->SetClusMaxThresh(0); + clusterizer->SetClusSigThresh(0); + //clusterizer->SetVerbose(1); + clusterizer->Dump(); + + // Open the AliRoot file + clusterizer->Open(infile,0); + //clusterizer->Open(infile,outfile,0); + + // Load the digits + clusterizer->ReadDigits(); + + // Find the cluster + clusterizer->MakeClusters(); + + // Write the cluster tree into file AliTRDclusters.root + clusterizer->WriteClusters(-1); + + // Save the clusterizer class in the AliROOT file + // clusterizer->Write(); + +} diff --git a/TRD/AliTRDdisplayDigits3D.C b/TRD/AliTRDdisplayDigits3D.C index c60c1031b85..2a28d3a162c 100644 --- a/TRD/AliTRDdisplayDigits3D.C +++ b/TRD/AliTRDdisplayDigits3D.C @@ -1,5 +1,5 @@ //_____________________________________________________________________________ -Int_t AliTRDdisplayDigits3D(Int_t event = 0, Int_t thresh = 0 +Int_t AliTRDdisplayDigits3D(Int_t event = 0, Int_t thresh = 2 , Bool_t sdigits = kFALSE) { // diff --git a/TRD/AliTRDdisplayTracks.C b/TRD/AliTRDdisplayTracks.C new file mode 100644 index 00000000000..d38f1de943d --- /dev/null +++ b/TRD/AliTRDdisplayTracks.C @@ -0,0 +1,143 @@ +void AliTRDdisplayTracks(Int_t track) +{ + + c1 = new TCanvas( "RecPoints", "RecPoints Display", 10, 10, 710, 740); + c1->SetFillColor(1); + TView *v=new TView(1); + + v->SetRange(-350.,-350.,-400.,350.,350.,400.); // full + // v->SetRange(0.,0.,0.,350.,350.,400.); // top right + // v->SetRange(-350.,0.,0.,0.,350.,400.); // top left + // v->SetRange(0.,-350.,0.,350.,0.,400.); // bottom right + // v->SetRange(-350.,-350.,0.,0.,0.,400.); // bottom left + // v->Side(); + v->Top(); + cerr<<"psi = "<GetPsi()<GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + cout << "Loaded shared libraries" << endl; + } + + +// Load tracks + + TFile *tf=TFile::Open("AliTRDtracks.root"); + if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return 3;} + TObjArray tarray(2000); + TTree *tracktree=(TTree*)tf->Get("TreeT"); + TBranch *tbranch=tracktree->GetBranch("tracks"); + Int_t nentr=tracktree->GetEntries(); + cerr<<"Found "<SetAddress(&iotrack); + tracktree->GetEvent(i); + tarray.AddLast(iotrack); + cerr<<"got track "<ReadClusters(RecPointsArray,alifile,-1); + Int_t nRecPoints = RecPointsArray->GetEntriesFast(); + cerr<<"Found "<GetListOfFiles()->FindObject(alifile); + if (!gafl) { + cout << "Open the ALIROOT-file " << alifile << endl; + gafl = new TFile(alifile); + } + else { + cout << alifile << " is already open" << endl; + } + + // Get AliRun object from file or create it if not on file + gAlice = (AliRun*) gafl->Get("gAlice"); + if (gAlice) + cout << "AliRun object found on file" << endl; + else + gAlice = new AliRun("gAlice","Alice test program"); + + + AliTRDv1 *TRD = (AliTRDv1*) gAlice->GetDetector("TRD"); + AliTRDgeometry *fGeom = TRD->GetGeometry(); + + Int_t i,j,index,det,sector, ti[3]; + Double_t x,y,z, cs,sn,tmp; + Float_t global[3], local[3]; + + // Display all TRD RecPoints + TPolyMarker3D *pm = new TPolyMarker3D(nRecPoints); + + for (Int_t i = 0; i < nRecPoints; i++) { + printf("\r point %d out of %d",i,nRecPoints); + AliTRDrecPoint *rp = (AliTRDrecPoint *) RecPointsArray->UncheckedAt(i); + + local[0]=rp->GetLocalRow(); + local[1]=rp->GetLocalCol(); + local[2]=rp->GetLocalTime(); + det=rp->GetDetector(); + + for (j = 0; j < 3; j++) { ti[j] = rp->GetTrackIndex(j); } + + if((track < 0) || + ((ti[0]==track)||(ti[1]==track)||(ti[2]==track))) { + + if(fGeom->Local2Global(det,local,global)) { + x=global[0]; y=global[1]; z=global[2]; + pm->SetPoint(i,x,y,z); + } + } + } + + pm->SetMarkerSize(1); pm->SetMarkerColor(10); pm->SetMarkerStyle(1); + pm->Draw(); + + + Int_t ntracks = tarray.GetEntriesFast(); + + for (i = 0; i < ntracks; i++) { + AliTRDtrack *pt = (AliTRDtrack *) tarray.UncheckedAt(i), &t=*pt; + Int_t nclusters = t.GetNclusters(); + cerr<<"in track "<UncheckedAt(index); + local[0]=rp->GetLocalRow(); + local[1]=rp->GetLocalCol(); + local[2]=rp->GetLocalTime(); + det=rp->GetDetector(); + + if(fGeom->Local2Global(det,local,global)) { + x=global[0]; y=global[1]; z=global[2]; + pm->SetPoint(j,x,y,z); + } + } + pm->SetMarkerSize(1); pm->SetMarkerColor(i%6+3); pm->SetMarkerStyle(1); + pm->Draw(); + } + + TGeometry *geom=(TGeometry*)gafl->Get("AliceGeom"); + geom->Draw("same"); + + c1->Modified(); + + c1->Update(); + + gafl->Close(); + +} diff --git a/TRD/AliTRDhits2digits.C b/TRD/AliTRDhits2digits.C new file mode 100644 index 00000000000..28a397000ba --- /dev/null +++ b/TRD/AliTRDhits2digits.C @@ -0,0 +1,47 @@ +void AliTRDhits2digits() +{ + +///////////////////////////////////////////////////////////////////////// +// +// Creates digits from the hit information. +// +///////////////////////////////////////////////////////////////////////// + + // Dynamically link some shared libs + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + cout << "Loaded shared libraries" << endl; + } + + // Input (and output) file name + Char_t *alifile = "galice.root"; + + // Create the TRD digitzer + AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer" + ,"TRD digitizer class"); + + // Set the parameter + digitizer->SetDebug(1); + + // Open the AliRoot file + digitizer->Open(alifile); + + + // 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 + digitizer->MakeDigits(); + + // Write the digits into the input file + digitizer->WriteDigits(); + + // Save the parameter object in the AliROOT file + parameter->Write(); + +} diff --git a/TRD/AliTRDhits2sdigits.C b/TRD/AliTRDhits2sdigits.C new file mode 100644 index 00000000000..d1f9f58639f --- /dev/null +++ b/TRD/AliTRDhits2sdigits.C @@ -0,0 +1,49 @@ +void AliTRDhits2sdigits() +{ + + ///////////////////////////////////////////////////////////////////////// + // + // Creates summable digits from the hit information. + // + ///////////////////////////////////////////////////////////////////////// + + // Dynamically link some shared libs + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + cout << "Loaded shared libraries" << endl; + } + + // Input (and output) file name + Char_t *alifile = "galice.root"; + + // Create the TRD digitzer + AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer" + ,"TRD digitizer class"); + + // Set the parameter + digitizer->SetDebug(1); + + // For the summable digits + digitizer->SetSDigits(kTRUE); + + // Open the AliRoot file + digitizer->Open(alifile); + + // 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 + digitizer->MakeDigits(); + + // Write the digits into the input file + digitizer->WriteDigits(); + + // Save the parameter object in the AliROOT file + parameter->Write(); + +} diff --git a/TRD/AliTRDmerge.C b/TRD/AliTRDmerge.C new file mode 100644 index 00000000000..2daa903009c --- /dev/null +++ b/TRD/AliTRDmerge.C @@ -0,0 +1,27 @@ +void AliTRDmerge() +{ + ///////////////////////////////////////////////////////////////////////// + // + // Merges different file with summable digits + // + ///////////////////////////////////////////////////////////////////////// + + AliRunDigitizer *manager = new AliRunDigitizer(2,1); + manager->SetInputStream(0,"galice_signal.root"); + manager->SetInputStream(1,"galice_background.root"); + + AliTRDdigitizer *digitizer = new AliTRDdigitizer(manager + ,"TRDdigitizer" + ,"TRD digitizer class"); + + // 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); + + // Do the merging + manager->Exec("deb"); + +} diff --git a/TRD/AliTRDrunSimple.C b/TRD/AliTRDrunSimple.C new file mode 100644 index 00000000000..6658c03d591 --- /dev/null +++ b/TRD/AliTRDrunSimple.C @@ -0,0 +1,274 @@ + +#if !defined(__CINT__) || defined(__MAKECINT__) + +#include "TH1.h" +#include "TFile.h" +#include "TCanvas.h" +#include "TProfile.h" +#include "TClonesArray.h" + +#include "AliModule.h" +#include "AliFRAMEv1.h" + +#include "AliTRDv1.h" +#include "AliTRDhit.h" +#include "AliTRDsim.h" +#include "AliTRDsimple.h" +#include "AliTRDsimpleGen.h" +#include "AliTRDsimpleMC.h" +#include "AliTRDgeometry.h" +#include "AliTRDdigitizer.h" +#include "AliTRDdigitsManager.h" +#include "AliTRDdataArrayI.h" + +#endif + +void AliTRDrunSimple() +{ + + /////////////////////////////////////////////////////////////////////////////// + // + // Macro to run the simplified simulator + // + /////////////////////////////////////////////////////////////////////////////// + + //_____________________________________________________________________________ + // + // Initialization part + //_____________________________________________________________________________ + // + + // The simple simulator + AliTRDsimple *simple = new AliTRDsimple(); + simple->Init(); + + // Initialize a dummy frame so that the TRD works + new AliFRAMEv1("FRAME","Space Frame"); + + // Initialize the TRD detector + AliTRDv1 *trd = new AliTRDv1("TRD","TRD slow simulator"); + + // Needed for some material properties + trd->CreateMaterials(); + + // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 85% Xe + 15% CO2) + trd->SetGasMix(1); + + // Get the pointer to the geometry object + AliTRDgeometry *geometry = trd->GetGeometry(); + + // The number of timebins + geometry->SetNTimeBin(30); + + // The additional timebins before and after the drift region + geometry->SetExpandTimeBin(5,15); + + // Switch on TR production + AliTRDsim *trdsim = trd->CreateTR(); + trdsim->SetNFoils(100); + trdsim->SetFoilThick(0.0013); + trdsim->SetGapThick(0.0500); + + // Initialize the TRD object + trd->Init(); + + // The digitizer + AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer","Digitizer class"); + digitizer->SetSimple(); + digitizer->SetGasGain(3200.); // GEANT + new Ermilova spectrum + digitizer->SetChipGain(6.1); + digitizer->SetNoise(1000.); + digitizer->SetADCinRange(1000.); + digitizer->SetADCoutRange(1023.); + digitizer->SetADCthreshold(0); + digitizer->InitDetector(); + //digitizer->SetTimeResponse(0); + //digitizer->SetVerbose(1); + + // The event generator + AliTRDsimpleGen *generator = simple->GetGenerator(); + generator->SetMomentum(3.0,3.0); + generator->SetPdg(11); // Electrons + //generator->SetPdg(211); // Pions + + //_____________________________________________________________________________ + // + // Histograms + //_____________________________________________________________________________ + // + + Int_t timeMax = geometry->GetTimeTotal(); + Int_t adcRange = ((Int_t) digitizer->GetADCoutRange()); + + TH1F *hQ = new TH1F("hQ" ,"Charge per hit (all)" ,100,0.0,1000.0); + TH1F *hQdedx = new TH1F("hQdedx","Charge per hit (dedx)",100,0.0,1000.0); + TH1F *hQtr = new TH1F("hQtr ","Charge per hit (tr)" ,100,0.0,1000.0); + + TProfile *hQX = new TProfile("hQX" ,"Charge per hit vs X (all)" ,35,0.0,3.5); + TProfile *hQXdedx = new TProfile("hQXdedx","Charge per hit vs X (dedx)",35,0.0,3.5); + TProfile *hQXtr = new TProfile("hQXtr ","Charge per hit vs X (tr)" ,35,0.0,3.5); + + TH1F *hNstep = new TH1F("hNstep","Number of steps / cm" , 151,-0.5, 150.5); + TH1F *hNel = new TH1F("hNel" ,"Number of electrons / cm",1001,-0.5,1000.5); + + TH1F *hAmp = new TH1F("hAmp","Amplitude of the digits" + ,adcRange+1,-0.5,((Float_t) adcRange)+0.5); + TProfile *hAmpTime = new TProfile("hAmpTime","Amplitude vs timebin" + ,timeMax,-0.5,((Float_t) timeMax)-0.5); + + //_____________________________________________________________________________ + // + // Event loop + //_____________________________________________________________________________ + // + + // Number of events (tracks) + Int_t nEvent = 10000; + + Float_t x0 = geometry->GetTime0(0) - AliTRDgeometry::DrThick(); + + TClonesArray *hitsArray = trd->Hits(); + + for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) { + + if (!(iEvent % 100) && (iEvent)) { + printf("Event no. %d\n",iEvent); + } + + // Generate the hits for one track + simple->ProcessEvent(iEvent); + + // Analyze the hits + Float_t nElectrons = 0.0; + for (Int_t iHit = 0; iHit < hitsArray->GetEntries(); iHit++) { + + AliTRDhit *hit = (AliTRDhit *) hitsArray->UncheckedAt(iHit); + Int_t charge = TMath::Abs(hit->GetCharge()); + Float_t x = hit->X() - x0; + hQ->Fill(charge); + hQX->Fill(x,charge); + if (hit->FromDrift() || + hit->FromAmplification()) { + hQdedx->Fill(charge); + hQXdedx->Fill(x,charge); + nElectrons += charge; + } + if (hit->FromTRphoton()) { + hQtr->Fill(charge); + hQXtr->Fill(x,charge); + } + + } + + hNstep->Fill(((AliTRDsimpleMC *) gMC)->GetNStep() / 3.5); + hNel->Fill(nElectrons / 3.5); + + // Digitize the hits + digitizer->MakeDigits(); + + // Analyze the digits + Int_t row = 8; + Int_t col = 48; + AliTRDdigitsManager *digitsManager = digitizer->Digits(); + AliTRDdataArrayI *digitsArray = digitsManager->GetDigits(12); + for (Int_t time = 0; time < timeMax; time++) { + + Int_t amp = digitsArray->GetDataUnchecked(row,col,time); + if (amp > 0) { + hAmp->Fill((Double_t) amp); + hAmpTime->Fill((Double_t) time, (Double_t) amp); + } + + } + + } + + //_____________________________________________________________________________ + // + // Plot the spectra + //_____________________________________________________________________________ + // + + TCanvas *cHit = new TCanvas("cHit","Hits",50,50,800,600); + cHit->Divide(3,2); + + cHit->cd(1); + gPad->SetLogy(); + hQ->SetLineColor(4); + hQ->SetXTitle("Q (a.u.)"); + hQ->SetYTitle("Entries"); + hQ->Draw(); + hQdedx->SetLineColor(3); + hQdedx->Draw("SAME"); + hQtr->SetLineColor(2); + hQtr->Draw("SAME"); + + cHit->cd(2); + gPad->SetLogy(); + hQX->SetLineColor(4); + hQX->SetXTitle("x (cm)"); + hQX->SetYTitle(" (a.u.)"); + hQX->Draw("HIST"); + hQXdedx->SetLineColor(3); + hQXdedx->Draw("SAMEHIST"); + + cHit->cd(3); + gPad->SetLogy(); + hQXtr->SetLineColor(2); + hQXtr->SetXTitle("x (cm)"); + hQXtr->SetYTitle(" (a.u.)"); + hQXtr->Draw("HIST"); + + cHit->cd(4); + hNstep->SetLineColor(4); + hNstep->SetXTitle("N_{step} / cm"); + hNstep->SetYTitle("Entries"); + hNstep->Draw(); + + cHit->cd(5); + hNel->SetLineColor(4); + hNel->SetXTitle("N_{electron} / cm"); + hNel->SetYTitle("Entries"); + hNel->Draw(); + + TCanvas *cDigit = new TCanvas("cDigit","Digits",50,50,600,400); + cDigit->Divide(2,1); + + cDigit->cd(1); + gPad->SetLogy(); + hAmp->SetLineColor(2); + hAmp->SetXTitle("ADC channels"); + hAmp->SetYTitle("Entries"); + hAmp->Draw(); + + cDigit->cd(2); + hAmpTime->SetLineColor(2); + hAmpTime->SetXTitle("Timebin"); + hAmpTime->SetYTitle(""); + hAmpTime->Draw("HIST"); + + //_____________________________________________________________________________ + // + // Save the histograms + //_____________________________________________________________________________ + // + + TFile *fileOut = new TFile("simple.root","RECREATE"); + + hQ->Write(); + hQdedx->Write(); + hQtr->Write(); + + hQX->Write(); + hQXdedx->Write(); + hQXtr->Write(); + + hNstep->Write(); + hNel->Write(); + + hAmp->Write(); + hAmpTime->Write(); + + fileOut->Close(); + +} diff --git a/TRD/AliTRDsdigits2digits.C b/TRD/AliTRDsdigits2digits.C new file mode 100644 index 00000000000..4f82ac90840 --- /dev/null +++ b/TRD/AliTRDsdigits2digits.C @@ -0,0 +1,54 @@ +void AliTRDsdigits2digits() +{ + + ///////////////////////////////////////////////////////////////////////// + // + // Converts s-digits to normal digits + // + ///////////////////////////////////////////////////////////////////////// + + // Dynamically link some shared libs + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + cout << "Loaded shared libraries" << endl; + } + + Char_t *fileName = "galice.root"; + + // Create the TRD digits merger + AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer" + ,"TRD digitizer class"); + + // Set the parameter + digitizer->SetDebug(1); + + // Initialize the geometry + digitizer->Open(fileName); + + // 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 manager for the input s-digits + AliTRDdigitsManager *sdigitsManager = new AliTRDdigitsManager(); + sdigitsManager->SetDebug(1); + sdigitsManager->SetSDigits(kTRUE); + sdigitsManager->ReadDigits(); + + // Add the s-digits to the input list + digitizer->AddSDigitsManager(sdigitsManager); + + // Convert the s-digits to normal digits + digitizer->SDigits2Digits(); + + // Store the digits + digitizer->WriteDigits(); + + // Save the parameter object in the AliROOT file + parameter->Write(); + +} diff --git a/TRD/AliTRDtest.C b/TRD/AliTRDtest.C index d7a07b283a5..cff2b548285 100644 --- a/TRD/AliTRDtest.C +++ b/TRD/AliTRDtest.C @@ -18,11 +18,9 @@ Int_t AliTRDtest() gAlice = (AliRun *) file->Get("gAlice"); // Analyze the TRD hits - gROOT->LoadMacro("$(ALICE_ROOT)/TRD/AliTRDanalyzeHits.C"); if (rc = AliTRDanalyzeHits()) return rc; // Run the digitization - gROOT->LoadMacro("$(ALICE_ROOT)/TRD/AliTRDcreateDigits.C"); if (rc = AliTRDcreateDigits()) return rc; if (gAlice) delete gAlice; @@ -30,11 +28,9 @@ Int_t AliTRDtest() gAlice = (AliRun *) file->Get("gAlice"); // Analyze the digits - gROOT->LoadMacro("$(ALICE_ROOT)/TRD/AliTRDanalyzeDigits.C"); if (rc = AliTRDanalyzeDigits()) return rc; // Create the cluster - gROOT->LoadMacro("$(ALICE_ROOT)/TRD/AliTRDcreateCluster.C"); if (rc = AliTRDcreateCluster()) return rc; if (gAlice) delete gAlice; @@ -42,9 +38,686 @@ Int_t AliTRDtest() gAlice = (AliRun *) file->Get("gAlice"); // Analyze the digits - gROOT->LoadMacro("$(ALICE_ROOT)/TRD/AliTRDanalyzeCluster.C"); if (rc = AliTRDanalyzeCluster()) return rc; return rc; } + +//_____________________________________________________________________________ +Int_t AliTRDanalyzeHits() +{ + // + // Analyzes the hits and fills QA-histograms + // + + Int_t rc = 0; + + if (!gAlice) { + cout << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " 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 << " No tree with clusters found" << endl; + rc = 4; + return rc; + } + + // Get the pointer to the hit container + TObjArray *clusterArray = trd->RecPoints(); + if (!(clusterArray)) { + cout << " No clusterArray found" << endl; + rc = 5; + return rc; + } + + // Set the branch address + clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray); + Int_t nEntries = clusterTree->GetEntries(); + cout << " 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 << " Found " << countCluster << " cluster in total" << endl; + cout << " 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; + +} -- 2.39.3