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
18 //if (gAlice) delete gAlice;
19 //file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
20 //gAlice = (AliRun *) file->Get("gAlice");
22 // Analyze the TRD hits
23 if (rc = AliTRDanalyzeHits()) return rc;
25 // Run the digitization
26 if (rc = AliTRDcreateDigits()) return rc;
28 // if (gAlice) delete gAlice;
29 // file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
30 // gAlice = (AliRun *) file->Get("gAlice");
33 if (rc = AliTRDanalyzeDigits()) return rc;
35 // // Create the cluster
36 // if (rc = AliTRDcreateCluster()) return rc;
38 // if (gAlice) delete gAlice;
39 // file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
40 // gAlice = (AliRun *) file->Get("gAlice");
42 // // Analyze the cluster
43 // if (rc = AliTRDanalyzeCluster()) return rc;
45 //file = (TFile *) gROOT->GetListOfFiles()->FindObject("TRD_test.root");
49 //if (rc = AliTRDcreateTracks()) return rc;
55 //_____________________________________________________________________________
56 Int_t AliTRDanalyzeHits()
59 // Analyzes the hits and fills QA-histograms
64 AliRunLoader *rl = gAlice->GetRunLoader();
66 cout << "<AliTRDanalyzeHits> No RunLoader found" << endl;
71 // Import the Trees for the event nEvent in the file
74 AliLoader* loader = rl->GetLoader("TRDLoader");
76 cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl;
81 // Get the pointer to the TRD detector
82 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
84 cout << "<AliTRDanalyzeHits> No TRD detector found" << endl;
89 // Get the pointer to the geometry object
90 AliTRDgeometryFull *geo;
92 geo = (AliTRDgeometryFull *) trd->GetGeometry();
95 cout << "<AliTRDanalyzeHits> No TRD geometry found" << endl;
100 AliTRDCommonParam *par = AliTRDCommonParam::Instance();
102 // Define the histograms
103 TH1F *hQdedx = new TH1F("hQdedx","Charge dedx-hits",100,0.0,1000.0);
104 TH1F *hQtr = new TH1F("hQtr" ,"Charge TR-hits" ,100,0.0,1000.0);
106 Float_t rmin = geo->Rmin();
107 Float_t rmax = geo->Rmax();
108 Float_t length = geo->GetChamberLength(0,2);
109 Float_t width = geo->GetChamberWidth(0);
110 Int_t ncol = par->GetColMax(0);
111 Int_t nrow = par->GetRowMax(0,2,13);
112 Int_t ntime = ((Int_t) (rmax - rmin) / 22.0);
114 TH2F *hZY = new TH2F("hZY" ,"Y vs Z (chamber 0)", nrow,-length/2.,length/2.
116 TH2F *hXZ = new TH2F("hXZ" ,"Z vs X (plane 0)" , ncol, -width/2., width/2.
117 , nrow,-length/2.,length/2.);
119 // Get the pointer hit tree
120 TTree *hitTree = loader->TreeH();
122 cout << "<AliTRDanalyzeHits> No hit tree found" << endl;
130 // Get the number of entries in the hit tree
131 // (Number of primary particles creating a hit somewhere)
132 Int_t nTrack = (Int_t) hitTree->GetEntries();
133 cout << "<AliTRDanalyzeHits> Found " << nTrack
134 << " primary particles with hits" << endl;
136 // Loop through all entries in the tree
137 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
140 nBytes += hitTree->GetEvent(iTrack);
142 // Loop through the TRD hits
144 AliTRDhit *hit = (AliTRDhit *) trd->FirstHit(-1);
150 Float_t x = hit->X();
151 Float_t y = hit->Y();
152 Float_t z = hit->Z();
153 Float_t q = hit->GetCharge();
154 Int_t track = hit->Track();
155 Int_t det = hit->GetDetector();
156 Int_t plane = geo->GetPlane(det);
166 hQtr->Fill(TMath::Abs(q));
169 hit = (AliTRDhit *) trd->NextHit();
175 cout << "<AliTRDanalyzeHits> Found " << countHits << " hits in total" << endl;
177 TCanvas *cHits = new TCanvas("cHits","AliTRDanalyzeHits",50,50,600,600);
194 //_____________________________________________________________________________
195 Int_t AliTRDcreateDigits()
198 // Creates the digits from the hits of the slow simulator
204 cout << "<AliTRDcreateDigits> No AliRun object found" << endl;
210 // Create the TRD digitzer
211 AliTRDdigitizer *digitizer = new AliTRDdigitizer("TRDdigitizer"
212 ,"TRD digitizer class");
213 digitizer->Open("TRD_test.root");
214 digitizer->InitDetector();
215 digitizer->InitOutput(0);
218 digitizer->SetDebug(1);
220 // Define the parameter object
221 // If no external parameter object is defined,
222 // default parameter will be used
223 AliTRDparameter *parameter = new AliTRDparameter("TRDparameter"
224 ,"TRD parameter class");
225 digitizer->SetParameter(parameter);
228 if (!(digitizer->MakeDigits())) {
233 // Write the digits into the input file
234 //if (!(digitizer->MakeBranch())) {
239 // Write the digits into the input file
240 if (!(digitizer->WriteDigits())) {
245 // Save the parameter object in the AliROOT file
246 if (!(parameter->Write())) {
255 //_____________________________________________________________________________
256 Int_t AliTRDanalyzeDigits()
259 // Analyzes the digits
264 const Int_t kNpla = AliTRDgeometry::Nplan();
267 cout << "<AliTRDanalyzeDigits> No AliRun object found" << endl;
272 AliRunLoader *rl = gAlice->GetRunLoader();
274 cout << "<AliTRDanalyzeHits> No RunLoader found" << endl;
279 // Import the Trees for the event nEvent in the file
280 rl->LoadKinematics();
284 AliLoader* loader = rl->GetLoader("TRDLoader");
286 cout << "<AliTRDanalyzeHits> No TRDLoader found" << endl;
291 // Get the pointer to the TRD detector
292 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
294 cout << "<AliTRDanalyzeDigits> No TRD detector found" << endl;
299 // Get the parameter object
300 AliTRDSimParam *parSim = AliTRDSimParam::Instance();
301 AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
302 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
304 // Define the histograms
305 Int_t adcRange = ((Int_t) parSim->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 = cal->GetNumberOfTimeBins();
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 = parCom->GetRowMax(iPla,iCha,iSec);
362 Int_t colMax = parCom->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");