]>
Commit | Line | Data |
---|---|---|
409f8c84 | 1 | #include "iostream.h" |
2 | ||
3 | void SPDclusterTestDubna (Int_t evNumber1=0,Int_t evNumber2=0) | |
4 | { | |
5 | ///////////////////////////////////////////////////////////////////////// | |
6 | // This macro is a small example of a ROOT macro | |
7 | // illustrating how to read the output of GALICE | |
8 | // and do some analysis. | |
9 | // | |
10 | ///////////////////////////////////////////////////////////////////////// | |
11 | ||
12 | // Dynamically link some shared libs | |
13 | ||
14 | if (gClassTable->GetID("AliRun") < 0) { | |
15 | gROOT->LoadMacro("loadlibs.C"); | |
16 | loadlibs(); | |
17 | } else { | |
18 | delete gAlice; | |
19 | gAlice=0; | |
20 | } | |
21 | ||
22 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
23 | ||
24 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); | |
25 | if (!file) file = new TFile("galice.root"); | |
26 | file->ls(); | |
27 | ||
28 | // Get AliRun object from file or create it if not on file | |
29 | ||
30 | if (!gAlice) { | |
31 | gAlice = (AliRun*)file->Get("gAlice"); | |
32 | if (gAlice) printf("AliRun object found on file\n"); | |
33 | if (!gAlice) gAlice = new AliRun("gAlice","Alice test program"); | |
34 | } | |
35 | ||
36 | // ------------ Cluster and point analysis histogramms ------------ | |
37 | ||
38 | TH1F *Nxpix1 = new TH1F("Nxpix1","Cluster size in x(r*phi) direction for layer 1",20,0.,20.); | |
39 | TH1F *Nxpix2 = new TH1F("Nxpix2","Cluster size in x(r*phi) direction for layer 2",20,0.,20.); | |
40 | TH1F *Nzpix1 = new TH1F("Nzpix1","Cluster size in z direction for layer 1",15,0.,15.); | |
41 | TH1F *Nzpix2 = new TH1F("Nzpix2","Cluster size in z direction for layer 2",15,0.,15.); | |
42 | TH1F *Xpix1 = new TH1F("Xpix1","Local x coordinate (mm) for layer 1",20,-2.,18.); | |
43 | TH1F *Xpix2 = new TH1F("Xpix2","Local x coordinate (mm) for layer 2",20,-2.,18.); | |
44 | TH1F *Zpix1 = new TH1F("Zpix1","Local z coordinate (mm) for layer 1",90,-2.,88.); | |
45 | TH1F *Zpix2 = new TH1F("Zpix2","Lolac z coordinate (mm) for layer 2",90,-2.,88.); | |
46 | ||
47 | TH1F *Xres1 = new TH1F("Xres1","Xrec and Xgen difference (micr) for layers 1",100,-200.,200.); | |
48 | TH1F *Xres2 = new TH1F("Xres2","Xrec and Xgen difference (micr) for layers 2",100,-200.,200.); | |
49 | TH1F *Zres1 = new TH1F("Zres1","Zrec and Zgen difference (micr) for layers 1",100,-800.,800.); | |
50 | TH1F *Zres2 = new TH1F("Zres2","Zrec and Zgen difference (micr) for layers 2",100,-800.,800.); | |
51 | ||
52 | TH1F *Ptot1 = new TH1F("Ptot1","Total momentum (GeV/C) for layers 1",100,0.,5.); | |
53 | TH1F *Pz1 = new TH1F("Pz1","Pz (GeV/C) for layers 1",100,-5.,5.); | |
54 | TH1F *Theta1 = new TH1F("Theta1","Theta angle (rad) for layers 1",100,0.,4.); | |
55 | TH1F *Y1 = new TH1F("Y1","Rapidity for layers 1",100,-4.,4.); | |
56 | TH1F *Eta1 = new TH1F("Eta1","PseudoRapidity for layers 1",100,-4.,4.); | |
57 | TH1F *Y1Den = new TH1F("Y1Den","Rapidity for layers 1",100,-0.5,0.5); | |
58 | TH1F *Eta1Den = new TH1F("Eta1Den","PseudoRapidity for layers 1",100,-0.5,0.5); | |
59 | TH1F *Y1DenA = new TH1F("Y1DenA","Rapidity for layers 1",100,-0.5,0.5); | |
60 | TH1F *Eta1DenA = new TH1F("Eta1DenA","PseudoRapidity for layers 1",100,-0.5,0.5); | |
61 | TH1F *Phi1 = new TH1F("Phi1","Phi angle (rad) for layers 1",100,0.,7.); | |
62 | TH1F *Ptot2 = new TH1F("Ptot2","Total momentum (GeV/C) for layers 2",100,0.,5.); | |
63 | TH1F *Pz2 = new TH1F("Pz2","Pz (GeV/C) for layers 2",100,-5.,5.); | |
64 | TH1F *Theta2 = new TH1F("Theta2","Theta angle (rad) for layers 2",100,0.,4.); | |
65 | TH1F *Y2 = new TH1F("Y2","Rapidity for layers 2",100,-4.,4.); | |
66 | TH1F *Eta2 = new TH1F("Eta2","PseudoRapidity for layers 2",100,-4.,4.); | |
67 | TH1F *Y2Den = new TH1F("Y2Den","Rapidity for layers 2",100,-0.5,0.5); | |
68 | TH1F *Eta2Den = new TH1F("Eta2Den","PseudoRapidity for layers 2",100,-0.5,0.5); | |
69 | TH1F *Y2DenA = new TH1F("Y2DenA","Rapidity for layers 2",100,-0.5,0.5); | |
70 | TH1F *Eta2DenA = new TH1F("Eta2DenA","PseudoRapidity for layers 2",100,-0.5,0.5); | |
71 | TH1F *Phi2 = new TH1F("Phi2","Phi angle (rad) for layers 2",100,0.,7.); | |
72 | ||
73 | // -------------- Create ntuples -------------------- | |
74 | // ntuple structures: | |
75 | ||
76 | struct { | |
77 | Int_t lay; | |
78 | Int_t nx; | |
79 | Int_t nz; | |
80 | Int_t hitprim; | |
81 | Int_t partcode; | |
82 | Int_t ntrover; | |
83 | Float_t dx; | |
84 | Float_t dz; | |
85 | Float_t pmod; | |
86 | } ntuple_st; | |
87 | ||
88 | struct { | |
89 | Int_t lay; | |
90 | Int_t lad; | |
91 | Int_t det; | |
92 | Int_t nx; | |
93 | Int_t nz; | |
94 | Int_t ntrover; | |
95 | Int_t noverlaps; | |
96 | Int_t noverprim; | |
97 | Float_t qcl; | |
98 | Float_t dx; | |
99 | Float_t dz; | |
100 | Float_t x; | |
101 | Float_t z; | |
102 | } ntuple1_st; | |
103 | ||
104 | struct { | |
105 | // Int_t lay; | |
106 | Int_t lay; | |
107 | Int_t nx; | |
108 | Int_t nz; | |
109 | Float_t x; | |
110 | Float_t z; | |
111 | Float_t qcl; | |
112 | } ntuple2_st; | |
113 | ||
114 | ntuple = new TTree("ntuple","Demo ntuple"); | |
115 | ntuple->Branch("lay",&ntuple_st.lay,"lay/I"); | |
116 | ntuple->Branch("nx",&ntuple_st.nx,"nx/I"); | |
117 | ntuple->Branch("nz",&ntuple_st.nz,"nz/I"); | |
118 | ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I"); | |
119 | ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I"); | |
120 | ntuple->Branch("ntrover",&ntuple_st.ntrover,"ntrover/I"); | |
121 | ntuple->Branch("dx",&ntuple_st.dx,"dx/F"); | |
122 | ntuple->Branch("dz",&ntuple_st.dz,"dz/F"); | |
123 | ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F"); | |
124 | ||
125 | ntuple1 = new TTree("ntuple1","Demo ntuple1"); | |
126 | ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I"); | |
127 | ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I"); | |
128 | ntuple1->Branch("det",&ntuple1_st.det,"det/I"); | |
129 | ntuple1->Branch("nx",&ntuple1_st.nx,"nx/I"); | |
130 | ntuple1->Branch("nz",&ntuple1_st.nz,"nz/I"); | |
131 | ntuple1->Branch("qcl",&ntuple1_st.qcl,"qcl/F"); | |
132 | ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I"); | |
133 | ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I"); | |
134 | ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I"); | |
135 | ntuple1->Branch("x",&ntuple1_st.x,"x/F"); | |
136 | ntuple1->Branch("z",&ntuple1_st.z,"z/F"); | |
137 | ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F"); | |
138 | ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F"); | |
139 | ||
140 | ||
141 | ntuple2 = new TTree("ntuple2","Demo ntuple2"); | |
142 | // ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I"); | |
143 | ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I"); | |
144 | ntuple2->Branch("x",&ntuple2_st.x,"x/F"); | |
145 | ntuple2->Branch("z",&ntuple2_st.z,"z/F"); | |
146 | ntuple2->Branch("nx",&ntuple2_st.nx,"nx/I"); | |
147 | ntuple2->Branch("nz",&ntuple2_st.nz,"nz/I"); | |
148 | ntuple2->Branch("qcl",&ntuple2_st.qcl,"qcl/F"); | |
149 | ||
150 | // ------------------------------------------------------------------------ | |
151 | // | |
152 | // Loop over events | |
153 | // | |
154 | for (int nev=0; nev<= evNumber2; nev++) { | |
155 | Int_t nparticles = gAlice->GetEvent(nev); | |
156 | cout << "nev " << nev <<endl; | |
157 | cout << "nparticles " << nparticles <<endl; | |
158 | if (nev < evNumber1) continue; | |
159 | if (nparticles <= 0) return; | |
160 | ||
161 | TTree *TH = gAlice->TreeH(); | |
162 | Int_t ntracks = TH->GetEntries(); | |
163 | cout<<"ntracks "<<ntracks<<endl; | |
164 | ||
165 | // Get pointers to Alice detectors and Digits containers | |
166 | AliITS *ITS = (AliITS*)gAlice->GetModule("ITS"); | |
167 | TClonesArray *Particles = gAlice->Particles(); | |
168 | ||
169 | if (ITS) { | |
170 | // fill modules with sorted by module hits | |
171 | Int_t nmodules; | |
172 | ITS->InitModules(-1,nmodules); | |
173 | // ITS->FillModules(nev,-1,evNumber2,nmodules," "," "); | |
174 | ITS->FillModules(nev,evNumber2,nmodules," "," "); | |
175 | //get pointer to modules array | |
176 | TObjArray *ITSmodules = ITS->GetModules(); | |
177 | AliITShit *itsHit; | |
178 | ||
179 | // get the Tree for clusters | |
180 | ITS->GetTreeC(nev); | |
181 | TTree *TC=ITS->TreeC(); | |
182 | Int_t nent=TC->GetEntries(); | |
183 | printf("Found %d entries in the tree (must be one per module per event!)\n",nent); | |
184 | Int_t lay, lad, det; | |
185 | AliITSgeom *geom = ITS->GetITSgeom(); | |
338e4f06 | 186 | |
187 | AliITSDetType *iDetType=ITS->DetType(0); | |
188 | AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel(); | |
189 | printf("SPD dimensions %f %f %f \n",seg0->Dx(),seg0->Dz(),seg0->Dy()); | |
190 | printf("SPD npixels %d %d \n",seg0->Npz(),seg0->Npx()); | |
191 | printf("SPD pitches %d %d \n",seg0->Dpz(0),seg0->Dpx(0)); | |
192 | Float_t SPDlength = seg0->Dz(); | |
193 | Float_t SPDwidth = seg0->Dx(); | |
194 | Float_t SPDthickness = seg0->Dy(); | |
195 | Float_t xpitch = seg0->Dpx(0); | |
196 | if(SPDlength > 80000) SPDthickness = 150; | |
197 | Float_t ylim; | |
198 | if(SPDthickness < 200) { | |
199 | ylim = SPDthickness/2 - 4; | |
200 | }else{ | |
201 | ylim = SPDthickness/2 - 10; | |
202 | } | |
203 | ||
409f8c84 | 204 | for (Int_t idettype=0;idettype<3;idettype++) { |
205 | ||
206 | TClonesArray *ITSclusters = ITS->ClustersAddress(idettype); | |
207 | //printf ("ITSclusters %p \n",ITSclusters); | |
208 | ||
209 | if (idettype != 0) continue; | |
210 | ||
211 | Float_t occup1 = 0; | |
212 | Float_t occup2 = 0; | |
213 | ||
214 | // Module loop | |
215 | for (Int_t mod=0; mod<nent; mod++) { | |
216 | AliITSmodule *itsModule = (AliITSmodule*)ITSmodules->At(mod); | |
217 | geom->GetModuleId(mod,lay,lad,det); | |
218 | ||
219 | Int_t nhits = itsModule->GetNhits(); | |
220 | //if(nhits) printf("module nhits %d %d\n",mod,nhits); | |
221 | if(!nhits) continue; | |
222 | ||
223 | ITS->ResetClusters(); | |
224 | TC->GetEvent(mod); | |
225 | Int_t nclust = ITSclusters->GetEntries(); | |
226 | if (!nclust) continue; | |
227 | ||
228 | // cluster/hit loops | |
229 | //cout<<"mod,lay,nclust,nhits ="<<mod<<","<<lay<<","<<nclust<<","<<nhits<<endl; | |
230 | for (Int_t clu=0;clu<nclust;clu++) { | |
231 | itsclu = (AliITSRawClusterSPD*)ITSclusters->UncheckedAt(clu); | |
232 | ||
233 | Int_t noverlaps = 0; | |
234 | Int_t noverprim = 0; | |
235 | ||
236 | Int_t clustersizex = itsclu->NclX(); | |
237 | Int_t clustersizez = itsclu->NclZ(); | |
409f8c84 | 238 | Int_t xstart = itsclu->XStartf(); |
239 | Int_t xstop = itsclu->XStopf(); | |
338e4f06 | 240 | Float_t fxstart = xstart*xpitch; |
241 | Float_t fxstop = (xstop+1)*xpitch; | |
409f8c84 | 242 | Float_t zstart = itsclu->ZStart(); |
243 | Float_t zstop = itsclu->ZStop(); | |
244 | Int_t zend = itsclu->Zend(); | |
245 | Int_t ntrover = itsclu->NTracks(); | |
338e4f06 | 246 | Float_t clusterx = 0; |
409f8c84 | 247 | Float_t clusterz = itsclu->Z(); |
338e4f06 | 248 | |
249 | Int_t tr0, tr1, tr2; | |
250 | itsclu->GetTracks(tr0,tr1,tr2); | |
251 | ||
252 | for(Int_t ii=0;ii<clustersizex;ii++) { | |
253 | clusterx += (xstart+0.5+ii)*xpitch; | |
254 | } | |
255 | clusterx /= clustersizex; | |
256 | clusterz /= clustersizez; | |
257 | ||
409f8c84 | 258 | Float_t clusterQ = itsclu->Q(); |
259 | ||
260 | if(lay == 1) occup1 += clusterQ; | |
261 | if(lay == 2) occup2 += clusterQ; | |
262 | ||
263 | ntuple2_st.lay = lay; | |
264 | ntuple2_st.x = clusterx/1000.; | |
265 | ntuple2_st.z = clusterz/1000.; | |
266 | ntuple2_st.nx = clustersizex; | |
267 | ntuple2_st.nz = clustersizez; | |
268 | ntuple2_st.qcl = clusterQ; | |
269 | ||
270 | ntuple2->Fill(); | |
271 | ||
272 | Int_t icl = 0; | |
273 | Float_t dxprimlast = 10.e+6; | |
274 | Float_t dzprimlast = 10.e+6; | |
275 | ||
338e4f06 | 276 | Float_t xhit0 = 1e+6; |
277 | Float_t zhit0 = 1e+6; | |
409f8c84 | 278 | for (Int_t hit=0;hit<nhits;hit++) { |
279 | ||
280 | // Find coordinate differences between the hit and cluster positions | |
281 | // for the resolution determination. | |
282 | ||
283 | itsHit = (AliITShit*)itsModule->GetHit(hit); | |
284 | ||
285 | Int_t hitlayer = itsHit->GetLayer(); | |
286 | Int_t hitladder= itsHit->GetLadder(); | |
287 | Int_t hitdet= itsHit->GetDetector(); | |
288 | Float_t dEn = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV | |
289 | Int_t track = itsHit->GetTrack(); | |
290 | Int_t dray = 0; | |
291 | Int_t hitstat = itsHit->GetTrackStatus(); | |
409f8c84 | 292 | Float_t zhit = 10000*itsHit->GetZL(); |
293 | Float_t xhit = 10000*itsHit->GetXL(); | |
338e4f06 | 294 | Float_t yhit = 10000*itsHit->GetYL(); |
295 | Int_t parent = itsHit->GetParticle()->GetFirstMother(); | |
296 | Int_t partcode = itsHit->GetParticle()->GetPdgCode(); | |
297 | ||
298 | Int_t hitprim = 0; | |
299 | ||
300 | if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery | |
301 | // particles | |
409f8c84 | 302 | |
303 | Float_t pxsimL = itsHit->GetPXL(); // the momenta at GEANT points | |
304 | Float_t pysimL = itsHit->GetPYL(); | |
305 | Float_t pzsimL = itsHit->GetPZL(); | |
306 | Float_t psimL = TMath::Sqrt(pxsimL*pxsimL+pysimL*pysimL+pzsimL*pzsimL); | |
307 | ||
308 | // Check boundaries | |
309 | if(zhit > SPDlength/2) { | |
310 | //cout<<"!!! z outside ="<<zhit<<endl; | |
311 | zhit = SPDlength/2 - 10; | |
312 | } | |
313 | if(zhit < 0 && zhit < -SPDlength/2) { | |
314 | //cout<<"!!! z outside ="<<zhit<<endl; | |
315 | zhit = -SPDlength/2 + 10; | |
316 | } | |
317 | if(xhit > SPDwidth/2) { | |
318 | //cout<<"!!! x outside ="<<xhit<<endl; | |
319 | xhit = SPDwidth/2 - 10; | |
320 | } | |
321 | if(xhit < 0 && xhit < -SPDwidth/2) { | |
322 | //cout<<"!!! x outside ="<<xhit<<endl; | |
323 | xhit = -SPDwidth/2 + 10; | |
324 | } | |
325 | ||
326 | zhit += SPDlength/2; | |
327 | xhit += SPDwidth/2; | |
409f8c84 | 328 | |
338e4f06 | 329 | |
330 | ||
331 | if(hitlayer == 1 && hitstat == 66 && yhit > ylim) { | |
332 | //if(hitstat == 66 && hitprim ==1) { | |
409f8c84 | 333 | xhit0 = xhit; |
334 | zhit0 = zhit; | |
335 | } | |
338e4f06 | 336 | |
337 | ||
338 | if(hitlayer == 2 && hitstat == 66 && yhit < -ylim) { | |
409f8c84 | 339 | xhit0 = xhit; |
340 | zhit0 = zhit; | |
338e4f06 | 341 | } |
342 | ||
409f8c84 | 343 | if(hitstat != 68) continue; // Take only the hit if the last |
344 | // track point went out from | |
345 | // the detector. | |
346 | ||
338e4f06 | 347 | if(xhit0 > 1e+5 || zhit0 > 1e+5) continue; |
409f8c84 | 348 | |
349 | Float_t xmed = (xhit + xhit0)/2; | |
350 | Float_t zmed = (zhit + zhit0)/2; | |
351 | ||
352 | Float_t xdif = xmed - clusterx; | |
353 | Float_t zdif = zmed - clusterz; | |
338e4f06 | 354 | xhit0 = 1e+6; |
355 | zhit0 = 1e+6; | |
356 | ||
357 | // Consider the hits inside of cluster region only b.b. | |
409f8c84 | 358 | |
338e4f06 | 359 | if((xmed >= fxstart && xmed <= fxstop) && (zmed >= zstart && zmed <= zstop)) { |
409f8c84 | 360 | |
338e4f06 | 361 | // Consider the hits only with the track number equaled to one |
362 | // of the cluster | |
363 | //if((track == tr0) || (track == tr1) || (track == tr2)) { | |
409f8c84 | 364 | |
409f8c84 | 365 | icl = 1; |
366 | ||
367 | // part = (TParticle *)particles.UncheckedAt(track); | |
368 | // Int_t partcode = part->GetPdgCode(); | |
369 | // Int_t primery = gAlice->GetPrimary(track); | |
370 | ||
409f8c84 | 371 | |
372 | // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+ | |
373 | // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda | |
374 | ||
375 | Float_t pmod = itsHit->GetParticle()->P(); // total momentum at the | |
376 | // vertex | |
377 | Float_t energy = itsHit->GetParticle()->Energy(); // energy at the | |
378 | // vertex | |
379 | Float_t mass = itsHit->GetParticle()->GetMass(); // particle mass | |
380 | ||
381 | Float_t pz = itsHit->GetParticle()->Pz(); // z momentum componetnt | |
382 | // at the vertex | |
383 | Float_t px = itsHit->GetParticle()->Px(); // z momentum componetnt | |
384 | // at the vertex | |
385 | Float_t py = itsHit->GetParticle()->Py(); // z momentum componetnt | |
386 | // at the vertex | |
387 | Float_t phi = itsHit->GetParticle()->Phi(); // Phi angle at the | |
388 | // vertex | |
389 | Float_t theta = itsHit->GetParticle()->Theta(); // Theta angle at the | |
390 | // vertex | |
391 | //Float_t eta = itsHit->GetParticle()->Eta(); // Pseudo rapidity at the | |
338e4f06 | 392 | // vertex |
393 | Float_t y; | |
409f8c84 | 394 | if((energy-pz) > 0) { |
338e4f06 | 395 | y = 0.5*TMath::Log((energy+pz)/(energy-pz)); |
409f8c84 | 396 | }else{ |
397 | cout<<" Warning: energy < pz ="<<energy<<","<<pz<<endl; | |
398 | y = 10; | |
399 | } | |
400 | Float_t eta = -TMath::Log(TMath::Tan(theta/2)); | |
401 | pmod *= 1.0e+3; | |
402 | ||
403 | ||
404 | Float_t pxsim = itsHit->GetPXG(); // the momenta at this GEANT point | |
405 | Float_t pysim = itsHit->GetPYG(); | |
406 | Float_t pzsim = itsHit->GetPZG(); | |
407 | Float_t psim = TMath::Sqrt(pxsim*pxsim+pysim*pysim+pzsim*pzsim); | |
408 | ||
338e4f06 | 409 | //Int_t hitprim = 0; |
409f8c84 | 410 | |
411 | if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e- | |
412 | // at p < 6 MeV/c | |
413 | ||
414 | if(dray == 0) noverlaps = noverlaps + 1; // overlapps for all hits but | |
415 | // not for delta ray which | |
416 | // also went out from the | |
417 | // detector and returned | |
418 | // again | |
419 | ||
338e4f06 | 420 | //if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery |
409f8c84 | 421 | // particles |
422 | ||
423 | if(hitprim > 0) noverprim = noverprim + 1; | |
424 | ||
425 | if(hitprim > 0) { | |
426 | dxprimlast = xdif; | |
427 | dzprimlast = zdif; | |
428 | } | |
429 | ||
430 | // fill ntuple | |
431 | ||
432 | ntuple_st.lay = hitlayer; | |
433 | ntuple_st.nx = clustersizex; | |
434 | ntuple_st.nz = clustersizez; | |
435 | ntuple_st.hitprim = hitprim; | |
436 | ntuple_st.partcode = partcode; | |
437 | ntuple_st.ntrover = ntrover; | |
438 | ntuple_st.dx = xdif; | |
439 | ntuple_st.dz = zdif; | |
440 | ntuple_st.pmod = pmod; | |
441 | ||
442 | ntuple->Fill(); | |
443 | ||
444 | if(hitlayer == 1) { | |
445 | Y1DenA->Fill(y); | |
446 | Eta1DenA->Fill(eta); | |
447 | } | |
448 | if(hitlayer == 2) { | |
449 | Y2DenA->Fill(y); | |
450 | Eta2DenA->Fill(eta); | |
451 | } | |
452 | ||
453 | ||
454 | ||
455 | if(hitprim > 0) { // for primary particles | |
456 | ||
457 | if(hitlayer == 1) { | |
458 | Xres1->Fill(xdif); | |
459 | Zres1->Fill(zdif); | |
460 | Ptot1->Fill(pmod/1000.); | |
461 | Pz1->Fill(pz); | |
462 | Theta1->Fill(theta); | |
463 | Y1->Fill(y); | |
464 | Eta1->Fill(eta); | |
465 | Y1Den->Fill(y); | |
466 | Eta1Den->Fill(eta); | |
467 | Phi1->Fill(phi); | |
468 | } | |
469 | if(hitlayer == 2) { | |
470 | Xres2->Fill(xdif); | |
471 | Zres2->Fill(zdif); | |
472 | Ptot2->Fill(pmod/1000.); | |
473 | Pz2->Fill(pz); | |
474 | Theta2->Fill(theta); | |
475 | Y2->Fill(y); | |
476 | Eta2->Fill(eta); | |
477 | Y2Den->Fill(y); | |
478 | Eta2Den->Fill(eta); | |
479 | Phi2->Fill(phi); | |
480 | } | |
481 | } // primery particles | |
482 | ||
483 | } // end of cluster region | |
338e4f06 | 484 | |
409f8c84 | 485 | } // end of hit loop |
486 | ||
487 | if(icl == 1) { | |
488 | ||
489 | // fill ntuple1 | |
490 | ||
409f8c84 | 491 | if(noverlaps == 0) noverlaps = 1; // cluster contains one or more |
492 | // delta rays only | |
493 | ||
494 | ntuple1_st.lay = lay; | |
495 | ntuple1_st.lad = lad; | |
496 | ntuple1_st.det = det; | |
497 | ntuple1_st.x = clusterx*1000.; | |
498 | ntuple1_st.z = clusterz*1000.; | |
499 | ntuple1_st.nx = clustersizex; | |
500 | ntuple1_st.nz = clustersizez; | |
501 | ntuple1_st.qcl = clusterQ; | |
502 | ntuple1_st.ntrover = ntrover; | |
503 | ntuple1_st.noverlaps = noverlaps; | |
504 | ntuple1_st.noverprim = noverprim; | |
505 | ntuple1_st.dx = dxprimlast; | |
506 | ntuple1_st.dz = dzprimlast; | |
507 | ||
508 | ntuple1->Fill(); | |
509 | ||
510 | } // icl = 1 | |
511 | } // cluster loop | |
512 | } // module loop | |
513 | ||
514 | cout<<" Occupancy for layer-1 ="<<occup1<<endl; | |
515 | cout<<" Occupancy for layer-2 ="<<occup2<<endl; | |
516 | // The real occupancy values are: | |
517 | // (for full ALICE event at the full SPD acceptence) | |
518 | // occup1 /= 3932160; | |
519 | // occup2 /= 7864320; | |
520 | ||
521 | } // idettype loop | |
522 | } // end if ITS | |
523 | } // event loop | |
524 | ||
525 | ||
526 | // Write and Draw Histogramms and ntuples | |
527 | ||
528 | ||
529 | ||
530 | TFile fhistos("SPD_his_dubna.root","RECREATE"); | |
531 | ||
532 | ntuple->Write(); | |
533 | ntuple1->Write(); | |
534 | ntuple2->Write(); | |
535 | ||
536 | Nxpix1->Write(); | |
537 | Nzpix1->Write(); | |
538 | Nxpix2->Write(); | |
539 | Nzpix2->Write(); | |
540 | ||
541 | Xpix1->Write(); | |
542 | Zpix1->Write(); | |
543 | Xpix2->Write(); | |
544 | Zpix2->Write(); | |
545 | ||
546 | Xres1->Write(); | |
547 | Zres1->Write(); | |
548 | Xres2->Write(); | |
549 | Zres2->Write(); | |
550 | ||
551 | Ptot1->Write(); | |
552 | Pz1->Write(); | |
553 | Theta1->Write(); | |
554 | Y1->Write(); | |
555 | Eta1->Write(); | |
556 | Y1Den->Write(); | |
557 | Eta1Den->Write(); | |
558 | Y1DenA->Write(); | |
559 | Eta1DenA->Write(); | |
560 | Phi1->Write(); | |
561 | ||
562 | Ptot2->Write(); | |
563 | Pz2->Write(); | |
564 | Theta2->Write(); | |
565 | Y2->Write(); | |
566 | Eta2->Write(); | |
567 | Y2Den->Write(); | |
568 | Eta2Den->Write(); | |
569 | Y2DenA->Write(); | |
570 | Eta2DenA->Write(); | |
571 | Phi2->Write(); | |
572 | ||
573 | fhistos.Close(); | |
574 | cout<<"!!! Histogramms and ntuples were written"<<endl; | |
575 | ||
576 | TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700); | |
577 | c1->Divide(2,2); | |
578 | c1->cd(1); | |
579 | gPad->SetFillColor(33); | |
580 | Xres1->SetFillColor(42); | |
581 | Xres1->Draw(); | |
582 | c1->cd(2); | |
583 | gPad->SetFillColor(33); | |
584 | Zres1->SetFillColor(46); | |
585 | Zres1->Draw(); | |
586 | c1->cd(3); | |
587 | gPad->SetFillColor(33); | |
588 | Xres2->SetFillColor(42); | |
589 | Xres2->Draw(); | |
590 | c1->cd(4); | |
591 | gPad->SetFillColor(33); | |
592 | Zres2->SetFillColor(46); | |
593 | Zres2->Draw(); | |
594 | ||
595 | cout<<"END test for clusters and hits "<<endl; | |
596 | ||
597 | // file->Close(); | |
598 | } | |
599 | ||
600 | ||
601 |