--- /dev/null
+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();
+
+}
--- /dev/null
+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();
+
+}
--- /dev/null
+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();
+
+}
//_____________________________________________________________________________
-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)
{
//
--- /dev/null
+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 = "<<v->GetPsi()<<endl;
+
+ // Dynamically link some shared libs
+ if (gClassTable->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 "<<nentr<<" entries in the track tree"<<endl;
+ for (Int_t i=0; i<nentr; i++) {
+ AliTRDtrack *iotrack=new AliTRDtrack;
+ tbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(i);
+ tarray.AddLast(iotrack);
+ cerr<<"got track "<<i<<": index ="<<iotrack->GetLabel()<<endl;
+ }
+ tf->Close();
+
+
+ // Load clusters
+ Char_t *alifile = "AliTRDclusters.root";
+ Int_t nEvent = 0;
+ TObjArray rparray(2000);
+ TObjArray *RecPointsArray = &rparray;
+ AliTRDtracker *Tracker = new AliTRDtracker("dummy","dummy");
+ Tracker->ReadClusters(RecPointsArray,alifile,-1);
+ Int_t nRecPoints = RecPointsArray->GetEntriesFast();
+ cerr<<"Found "<<nRecPoints<<" rec. points"<<endl;
+
+
+ // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
+ alifile = "galice.root";
+ 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");
+
+
+ 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 "<<i<<" found "<<nclusters<<" clusters"<<endl;
+
+ TPolyMarker3D *pm = new TPolyMarker3D(nclusters);
+ for(j = 0; j < nclusters; j++) {
+ index = t.GetClusterIndex(j);
+ AliTRDrecPoint *rp = (AliTRDrecPoint *) RecPointsArray->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();
+
+}
--- /dev/null
+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();
+
+}
--- /dev/null
+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();
+
+}
--- /dev/null
+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");
+
+}
--- /dev/null
+
+#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("<Q> (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("<Q> (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("<ADC channels>");
+ 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();
+
+}
--- /dev/null
+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();
+
+}
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;
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;
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 << "<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;
+
+}