Changes in digits IO. Add merging of summable digits
[u/mrichter/AliRoot.git] / TRD / AliTRDanalyzeHits.C
CommitLineData
793ff80c 1Int_t AliTRDanalyzeHits()
2{
3 //
4 // Analyzes the hits and fills QA-histograms
5 //
6
7 Int_t rc = 0;
8
9 if (!gAlice) {
10 cout << "<AliTRDanalyzeHits> No AliRun object found" << endl;
11 rc = 1;
12 return rc;
13 }
14 gAlice->GetEvent(0);
15
16 // Get the pointer to the TRD detector
abaf1f1d 17 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
18 if (!trd) {
793ff80c 19 cout << "<AliTRDanalyzeHits> No TRD detector found" << endl;
20 rc = 2;
21 return rc;
22 }
23
f9428ca8 24 // Get the pointer to the geometry object
abaf1f1d 25 AliTRDgeometryFull *geo;
26 if (trd) {
27 geo = (AliTRDgeometryFull *) trd->GetGeometry();
f9428ca8 28 }
29 else {
30 cout << "<AliTRDanalyzeHits> No TRD geometry found" << endl;
31 rc = 3;
32 return rc;
33 }
34
793ff80c 35 // Define the histograms
f9428ca8 36 TH1F *hQdedx = new TH1F("hQdedx","Charge dedx-hits",100,0.0,1000.0);
37 TH1F *hQtr = new TH1F("hQtr" ,"Charge TR-hits" ,100,0.0,1000.0);
38
abaf1f1d 39 Float_t rmin = geo->Rmin();
40 Float_t rmax = geo->Rmax();
41 Float_t length = geo->GetChamberLengthI(0);
42 Float_t width = geo->GetChamberWidth(0);
43 Int_t ncol = geo->GetColMax(0);
44 Int_t nrow = geo->GetRowMax(0,2,13);
45 Int_t ntime = ((Int_t) (rmax - rmin) / geo->GetTimeBinSize());
f9428ca8 46
47 TH2F *hZY = new TH2F("hZY" ,"Y vs Z (chamber 0)", nrow,-length/2.,length/2.
48 ,ntime, rmin, rmax);
49 TH2F *hXZ = new TH2F("hXZ" ,"Z vs X (plane 0)" , ncol, -width/2., width/2.
50 , nrow,-length/2.,length/2.);
51
52 TH1F *hNprimP = new TH1F("hNprimP","Number of primary electrons per cm (MIP pion)"
53 ,50,0.0,100.0);
54 TH1F *hNprimE = new TH1F("hNprimE","Number of primary electrons per cm (3GeV electron)"
55 ,50,0.0,100.0);
56 TH1F *hNtotP = new TH1F("hNtotP" ,"Number of electrons per cm (MIP pion)"
57 ,50,0.0,1000.0);
58 TH1F *hNtotE = new TH1F("hNtotE" ,"Number of electrons per cm (3GeV electron)"
59 ,50,0.0,1000.0);
793ff80c 60
61 // Get the pointer hit tree
abaf1f1d 62 TTree *hitTree = gAlice->TreeH();
63 if (!hitTree) {
793ff80c 64 cout << "<AliTRDanalyzeHits> No hit tree found" << endl;
f9428ca8 65 rc = 4;
793ff80c 66 return rc;
67 }
68
69 Int_t countHits = 0;
70 Int_t nBytes = 0;
71
72 // Get the number of entries in the hit tree
73 // (Number of primary particles creating a hit somewhere)
abaf1f1d 74 Int_t nTrack = (Int_t) hitTree->GetEntries();
793ff80c 75 cout << "<AliTRDanalyzeHits> Found " << nTrack
76 << " primary particles with hits" << endl;
77
78 // Loop through all entries in the tree
79 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
80
81 gAlice->ResetHits();
abaf1f1d 82 nBytes += hitTree->GetEvent(iTrack);
793ff80c 83
f9428ca8 84 Int_t nPrimE = 0;
85 Int_t nPrimP = 0;
86 Int_t nTotE = 0;
87 Int_t nTotP = 0;
88
793ff80c 89 // Get the number of hits in the TRD created by this particle
abaf1f1d 90 Int_t nHit = trd->Hits()->GetEntriesFast();
793ff80c 91
92 // Loop through the TRD hits
93 for (Int_t iHit = 0; iHit < nHit; iHit++) {
94
95 countHits++;
96
abaf1f1d 97 AliTRDhit *hit = (AliTRDhit *) trd->Hits()->UncheckedAt(iHit);
793ff80c 98
f9428ca8 99 Float_t x = hit->X();
100 Float_t y = hit->Y();
101 Float_t z = hit->Z();
102 Float_t q = hit->GetCharge();
103 Int_t track = hit->Track();
104 Int_t det = hit->GetDetector();
abaf1f1d 105 Int_t plane = geo->GetPlane(det);
793ff80c 106
f9428ca8 107 if (hit->FromDrift()) {
793ff80c 108 hQdedx->Fill(q);
f9428ca8 109 hZY->Fill(z,y);
110 if (plane == 0) {
111 hXZ->Fill(x,z);
112 }
4a0fe73c 113 }
f9428ca8 114 else if (hit->FromTRphoton()) {
793ff80c 115 hQtr->Fill(TMath::Abs(q));
4a0fe73c 116 }
793ff80c 117
f9428ca8 118 TParticle *part = gAlice->Particle(track);
119
120 if ((plane == 0) && (hit->FromDrift())) {
121
122 // 3 GeV electrons
123 if ((part->GetPdgCode() == 11) &&
124 (part->P() > 2.9)) {
125 nPrimE++;
126 nTotE += ((Int_t) q);
127 }
128
129 // MIP pions
130 if ((part->GetPdgCode() == -211) &&
131 (part->P() > 0.5)) {
132 nPrimP++;
133 nTotP += ((Int_t) q);
134 }
135
136 }
793ff80c 137
138 }
139
f9428ca8 140 if (nPrimE > 0) hNprimE->Fill(((Double_t) nPrimE)/3.);
141 if (nPrimP > 0) hNprimP->Fill(((Double_t) nPrimP)/3.);
142 if (nTotE > 0) hNtotE->Fill(((Double_t) nTotE)/3.);
143 if (nTotP > 0) hNtotP->Fill(((Double_t) nTotP)/3.);
144
793ff80c 145 }
146
147 cout << "<AliTRDanalyzeHits> Found " << countHits << " hits in total" << endl;
148
149 TCanvas *cHits = new TCanvas("cHits","AliTRDanalyzeHits",50,50,600,600);
150 cHits->Divide(2,2);
151 cHits->cd(1);
152 hXZ->Draw("COL");
153 cHits->cd(2);
154 hZY->Draw("COL");
155 cHits->cd(3);
156 gPad->SetLogy();
157 hQdedx->Draw();
158 cHits->cd(4);
159 gPad->SetLogy();
160 hQtr->Draw();
161
f9428ca8 162 TCanvas *cNel = new TCanvas("cNel","Ionization",50,50,600,600);
163 cNel->Divide(2,2);
164 cNel->cd(1);
165 hNprimE->SetStats();
166 hNprimE->Draw();
167 cNel->cd(2);
168 hNprimP->SetStats();
169 hNprimP->Draw();
170 cNel->cd(3);
171 hNtotE->SetStats();
172 hNtotE->Draw();
173 cNel->cd(4);
174 hNtotP->SetStats();
175 hNtotP->Draw();
176
177 TFile *fout = new TFile("TRD_ionization.root","RECREATE");
178 hNprimE->Write();
179 hNprimP->Write();
180 hNtotE->Write();
181 hNtotP->Write();
182 fout->Close();
183
793ff80c 184 return rc;
185
186}