5 // Test macro for the TRD code
12 // Initialize the test setup
13 gAlice->Init("$(ALICE_ROOT)/TRD/AliTRDconfig.C");
15 // Run one event and create the hits
19 //if (gAlice) delete gAlice;
20 //file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
21 //gAlice = (AliRun *) file->Get("gAlice");
23 // Analyze the TRD hits
24 if (rc = AliTRDanalyzeHits()) return rc;
26 // Run the digitization
27 if (rc = AliTRDcreateDigits()) return rc;
29 // if (gAlice) delete gAlice;
30 // file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
31 // gAlice = (AliRun *) file->Get("gAlice");
34 if (rc = AliTRDanalyzeDigits()) return rc;
36 // // Create the cluster
37 // if (rc = AliTRDcreateCluster()) return rc;
39 // if (gAlice) delete gAlice;
40 // file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
41 // gAlice = (AliRun *) file->Get("gAlice");
43 // // Analyze the cluster
44 // if (rc = AliTRDanalyzeCluster()) return rc;
46 //file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
50 //if (rc = AliTRDcreateTracks()) return rc;
56 //_____________________________________________________________________________
57 Int_t AliTRDanalyzeHits()
60 // Analyzes the hits and fills QA-histograms
65 AliRunLoader *rl = gAlice->GetRunLoader();
67 cout << "<AliTRDanalyzeHits> No RunLoader found" << endl;
72 // Import the Trees for the event nEvent in the file
75 AliLoader* loader = rl->GetLoader("TRDLoader");
77 cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl;
82 // Get the pointer to the TRD detector
83 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
85 cout << "<AliTRDanalyzeHits> No TRD detector found" << endl;
90 // Get the pointer to the geometry object
91 AliTRDgeometryFull *geo;
93 geo = (AliTRDgeometryFull *) trd->GetGeometry();
96 cout << "<AliTRDanalyzeHits> No TRD geometry found" << endl;
101 AliTRDparameter *par = new AliTRDparameter("TRDparameter","TRD parameter class");
103 // Define the histograms
104 TH1F *hQdedx = new TH1F("hQdedx","Charge dedx-hits",100,0.0,1000.0);
105 TH1F *hQtr = new TH1F("hQtr" ,"Charge TR-hits" ,100,0.0,1000.0);
107 Float_t rmin = geo->Rmin();
108 Float_t rmax = geo->Rmax();
109 Float_t length = geo->GetChamberLength(0,2);
110 Float_t width = geo->GetChamberWidth(0);
111 Int_t ncol = par->GetColMax(0);
112 Int_t nrow = par->GetRowMax(0,2,13);
113 Int_t ntime = ((Int_t) (rmax - rmin) / par->GetTimeBinSize());
115 TH2F *hZY = new TH2F("hZY" ,"Y vs Z (chamber 0)", nrow,-length/2.,length/2.
117 TH2F *hXZ = new TH2F("hXZ" ,"Z vs X (plane 0)" , ncol, -width/2., width/2.
118 , nrow,-length/2.,length/2.);
120 // Get the pointer hit tree
121 TTree *hitTree = loader->TreeH();
123 cout << "<AliTRDanalyzeHits> No hit tree found" << endl;
131 // Get the number of entries in the hit tree
132 // (Number of primary particles creating a hit somewhere)
133 Int_t nTrack = (Int_t) hitTree->GetEntries();
134 cout << "<AliTRDanalyzeHits> Found " << nTrack
135 << " primary particles with hits" << endl;
137 // Loop through all entries in the tree
138 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
141 nBytes += hitTree->GetEvent(iTrack);
143 // Loop through the TRD hits
145 AliTRDhit *hit = (AliTRDhit *) trd->FirstHit(-1);
151 Float_t x = hit->X();
152 Float_t y = hit->Y();
153 Float_t z = hit->Z();
154 Float_t q = hit->GetCharge();
155 Int_t track = hit->Track();
156 Int_t det = hit->GetDetector();
157 Int_t plane = geo->GetPlane(det);
167 hQtr->Fill(TMath::Abs(q));
170 hit = (AliTRDhit *) trd->NextHit();
176 cout << "<AliTRDanalyzeHits> Found " << countHits << " hits in total" << endl;
178 TCanvas *cHits = new TCanvas("cHits","AliTRDanalyzeHits",50,50,600,600);
195 //_____________________________________________________________________________
196 Int_t AliTRDcreateDigits()
199 // Creates the digits from the hits of the slow simulator
205 cout << "<AliTRDcreateDigits> No AliRun object found" << endl;
211 // Create the TRD digitzer
212 AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
213 ,"TRD digitizer class");
214 digitizer->Open("TRD_test.root");
215 digitizer->InitDetector();
216 digitizer->InitOutput(0);
219 digitizer->SetDebug(1);
221 // Define the parameter object
222 // If no external parameter object is defined,
223 // default parameter will be used
224 AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
225 ,"TRD parameter class");
226 digitizer->SetParameter(parameter);
229 if (!(digitizer->MakeDigits())) {
234 // Write the digits into the input file
235 //if (!(digitizer->MakeBranch())) {
240 // Write the digits into the input file
241 if (!(digitizer->WriteDigits())) {
246 // Save the parameter object in the AliROOT file
247 if (!(parameter->Write())) {
256 //_____________________________________________________________________________
257 Int_t AliTRDanalyzeDigits()
260 // Analyzes the digits
265 const Int_t kNpla = AliTRDgeometry::Nplan();
268 cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
273 AliRunLoader *rl = gAlice->GetRunLoader();
275 cout << "<AliTRDanalyzeHits> No RunLoader found" << endl;
280 // Import the Trees for the event nEvent in the file
281 rl->LoadKinematics();
285 AliLoader* loader = rl->GetLoader("TRDLoader");
287 cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl;
292 // Get the pointer to the TRD detector
293 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
295 cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl;
300 // Get the parameter object
301 AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
302 ,"TRD parameter class");
304 // Define the histograms
305 Int_t adcRange = ((Int_t) parameter->GetADCoutRange());
306 TH1F *hAmpAll = new TH1F("hAmpAll" ,"Amplitude of the digits (all)"
307 ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
308 TH1F *hAmpEl = new TH1F("hAmpEl" ,"Amplitude of the digits (electrons)"
309 ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
310 TH1F *hAmpPi = new TH1F("hAmpPi" ,"Amplitude of the digits (pions)"
311 ,adcRange+1,-0.5,((Float_t) adcRange)+0.5);
312 TH1F *hAmpNoise = new TH1F("hAmpNoise","Amplitude of the digits (noise)"
315 // Get the pointer to the geometry object
318 geo = trd->GetGeometry();
321 cout << "<AliTRDanalyzeDigits> No TRD geometry found" << endl;
326 // Create the digits manager
327 AliTRDdigitsManager *digitsManager = new AliTRDdigitsManager();
328 digitsManager->SetDebug(1);
330 // Read the digits from the file
331 if (!(digitsManager->ReadDigits(loader->TreeD()))) {
332 cout << "<AliTRDanalyzeDigits> Cannot read the digits" << endl;
337 // Get the particle stack
338 AliStack *kineStack = rl->Stack();
340 cout << "<AliTRDanalyzeDigits> Cannot find the KINE stack" << endl;
345 AliTRDmatrix *matrix;
347 Int_t countDigits = 0;
348 Int_t iSec = trd->GetSensSector();
349 Int_t iCha = trd->GetSensChamber();
350 Int_t timeMax = parameter->GetTimeTotal();
352 TProfile *hAmpTimeEl = new TProfile("hAmpTimeEl","Amplitude of the digits (electrons)"
353 ,timeMax,-0.5,((Double_t) timeMax)-0.5);
354 TProfile *hAmpTimePi = new TProfile("hAmpTimePi","Amplitude of the digits (pions)"
355 ,timeMax,-0.5,((Double_t) timeMax)-0.5);
357 // Loop over all planes
358 for (Int_t iPla = 0; iPla < kNpla; iPla++) {
360 Int_t iDet = geo->GetDetector(iPla,iCha,iSec);
361 Int_t rowMax = parameter->GetRowMax(iPla,iCha,iSec);
362 Int_t colMax = parameter->GetColMax(iPla);
365 matrix = new AliTRDmatrix(rowMax,colMax,timeMax,iSec,iCha,iPla);
368 // Loop through the detector pixel
369 for (Int_t time = 0; time < timeMax; time++) {
370 for (Int_t col = 0; col < colMax; col++) {
371 for (Int_t row = 0; row < rowMax; row++) {
373 AliTRDdigit *digit = digitsManager->GetDigit(row,col,time,iDet);
374 Int_t amp = digit->GetAmp();
375 Int_t track0 = digitsManager->GetTrack(0,row,col,time,iDet);
376 Int_t track1 = digitsManager->GetTrack(1,row,col,time,iDet);
377 TParticle *particle = 0;
379 particle = (TParticle *) kineStack->Particle(track0);
385 matrix->SetSignal(row,col,time,amp);
394 hAmpNoise->Fill(amp);
398 if ((particle) && (particle->GetPdgCode() == 11) && (track1 < 0)) {
400 hAmpTimeEl->Fill(time,amp);
404 if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
406 hAmpTimePi->Fill(time,amp);
417 cout << "<AliTRDanalyzeDigits> Found " << countDigits << " digits in total" << endl;
419 // Display the detector matrix
422 TCanvas *cDigits = new TCanvas("cDigits","AliTRDanalyzeDigits",50,50,600,800);
423 cDigits->Divide(2,3);
426 hAmpAll->SetXTitle("Amplitude (ADC-channels)");
427 hAmpAll->SetYTitle("Entries");
431 hAmpNoise->SetXTitle("Amplitude (ADC-channels)");
432 hAmpNoise->SetYTitle("Entries");
436 hAmpEl->SetXTitle("Amplitude (ADC-channels)");
437 hAmpEl->SetYTitle("Entries");
441 hAmpPi->SetXTitle("Amplitude (ADC-channels)");
442 hAmpPi->SetYTitle("Entries");
445 hAmpTimeEl->SetXTitle("Timebin number");
446 hAmpTimeEl->SetYTitle("Mean amplitude");
447 hAmpTimeEl->Draw("HIST");
449 hAmpTimePi->SetXTitle("Timebin number");
450 hAmpTimePi->SetYTitle("Mean amplitude");
451 hAmpTimePi->Draw("HIST");
457 //_____________________________________________________________________________
458 Int_t AliTRDcreateCluster()
461 // Creates the cluster from the digits
466 // Create the clusterizer
467 AliTRDclusterizerV1 *clusterizer = new AliTRDclusterizerV1("TRDclusterizer"
468 ,"TRD clusterizer class");
469 clusterizer->SetVerbose(1);
471 // Define the parameter object
472 // If no external parameter object is defined,
473 // default parameter will be used
474 AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
475 ,"TRD parameter class");
476 parameter->SetClusMaxThresh(0);
477 parameter->SetClusSigThresh(0);
478 clusterizer->SetParameter(parameter);
481 if (!(clusterizer->Open("TRD_test.root",0))) {
487 if (!(clusterizer->ReadDigits())) {
493 if (!(clusterizer->MakeClusters())) {
498 // Write the cluster tree into the file
499 if (!(clusterizer->WriteClusters(-1))) {
504 // Save the clusterizer class in the file
505 if (!(clusterizer->Write())) {
514 //_____________________________________________________________________________
515 Int_t AliTRDanalyzeCluster()
518 // Analyzes the cluster
524 cout << "<AliTRDanalyzeCluster> No AliRun object found" << endl;
530 // Get the pointer to the TRD detector
531 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
533 cout << "<AliTRDanalyzeCluster> No TRD detector found" << endl;
538 // Define the histograms
539 TH1F *hClusAll = new TH1F("hClusAll" ,"Amplitude of the cluster (all)"
541 TH1F *hClusNoise = new TH1F("hClusNoise","Amplitude of the cluster (noise)"
543 TH1F *hClusEl = new TH1F("hClusEl" ,"Amplitude of the cluster (electron)"
545 TH1F *hClusPi = new TH1F("hClusPi" ,"Amplitude of the cluster (pion)"
548 // Get the pointer to the geometry object
551 geo = trd->GetGeometry();
554 cout << "<AliTRDanalyzeCluster> No TRD geometry found" << endl;
559 // Get the pointer to the hit-tree
560 TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
561 TTree *clusterTree = (TTree *) file->Get("TreeR0_TRD");
562 if (!(clusterTree)) {
563 cout << "<AliTRDanalyzeCluster> No tree with clusters found" << endl;
568 // Get the pointer to the hit container
569 TObjArray *clusterArray = trd->RecPoints();
570 if (!(clusterArray)) {
571 cout << "<AliTRDanalyzeCluster> No clusterArray found" << endl;
576 // Set the branch address
577 clusterTree->GetBranch("TRDcluster")->SetAddress(&clusterArray);
578 Int_t nEntries = clusterTree->GetEntries();
579 cout << "<AliTRDanalyzeCluster> Number of entries in the cluster tree = "
583 Int_t countCluster = 0;
584 Int_t countOverlap = 0;
586 // Loop through all entries in the tree
588 for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
591 nbytes += clusterTree->GetEvent(iEntry);
593 // Get the number of points in the detector
594 Int_t nCluster = clusterArray->GetEntriesFast();
596 // Loop through all TRD digits
597 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
599 // Get the information for this digit
600 AliTRDcluster *cluster = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
601 Int_t detector = cluster->GetDetector();
602 Int_t sector = geo->GetSector(detector);
603 Int_t plane = geo->GetPlane(detector);
604 Int_t chamber = geo->GetChamber(detector);
605 Float_t energy = cluster->GetQ();
606 Int_t track0 = cluster->GetLabel(0);
607 Int_t track1 = cluster->GetLabel(1);
608 Int_t track2 = cluster->GetLabel(2);
609 TParticle *particle = 0;
611 particle = gAlice->Particle(track0);
615 if (!cluster->Isolated()) countOverlap++;
618 hClusAll->Fill(energy);
620 if (cluster->Isolated()) {
624 hClusNoise->Fill(energy);
628 if ((particle) && (particle->GetPdgCode() == 11) && (track1 < 0)) {
629 hClusEl->Fill(energy);
633 if ((particle) && (particle->GetPdgCode() == -211) && (track1 < 0)) {
634 hClusPi->Fill(energy);
643 cout << "<AliTRDanalyzeCluster> Found " << countCluster << " cluster in total" << endl;
644 cout << "<AliTRDanalyzeCluster> Found " << countOverlap << " overlapping cluster" << endl;
647 TCanvas *cCluster = new TCanvas("cCluster","AliTRDanalyzeCluster",50,50,600,600);
648 cCluster->Divide(2,2);
653 hClusAll->Fit("landau","0");
654 fun = (TF1 *) hClusAll->GetListOfFunctions()->First();
655 Float_t meanAll = fun->GetParameter(1);
657 fun->SetLineColor(2);
662 Float_t meanNoise = hClusNoise->GetMean();
667 hClusEl->Fit("landau","0");
668 fun = (TF1 *) hClusEl->GetListOfFunctions()->First();
669 fun->SetLineColor(2);
670 Float_t meanEl = fun->GetParameter(1);
676 hClusPi->Fit("landau","0");
677 fun = (TF1 *) hClusPi->GetListOfFunctions()->First();
678 fun->SetLineColor(2);
679 Float_t meanPi = fun->GetParameter(1);
684 cout << "##################################################################" << endl;
685 cout << " Mean all = " << meanAll << endl;
686 cout << " Mean noise = " << meanNoise << endl;
687 cout << " Mean electrons = " << meanEl << endl;
688 cout << " Mean pions = " << meanPi << endl;
689 cout << "##################################################################" << endl;
696 //_____________________________________________________________________________
697 Int_t AliTRDcreateTracks()
700 // Creates the tracks
705 // Create the tracker
706 AliTRDtracker *tracker = new AliTRDtracker("TRDtracker","TRD tracker");
708 // Read in the kine tree and the cluster
709 tracker->GetEvent("TRD_test.root","TRD_test.root");
712 TH1F *hs = new TH1F("hs","hs",100,0.0,1.0);
713 TH1F *hd = new TH1F("hd","hd",100,0.0,1.0);
714 tracker->Clusters2Tracks(hs,hd);
717 tracker->WriteTracks("TRD_test.root");