3 void SPDclusterTestDubna (Int_t evNumber1=0,Int_t evNumber2=0)
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.
10 /////////////////////////////////////////////////////////////////////////
12 // Dynamically link some shared libs
14 if (gClassTable->GetID("AliRun") < 0) {
15 gROOT->LoadMacro("loadlibs.C");
22 // Connect the Root Galice file containing Geometry, Kine and Hits
24 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
25 if (!file) file = new TFile("galice.root");
28 // Get AliRun object from file or create it if not on file
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");
36 // ------------ Cluster and point analysis histogramms ------------
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.);
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.);
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.);
73 // -------------- Create ntuples --------------------
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");
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");
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");
150 // ------------------------------------------------------------------------
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;
161 TTree *TH = gAlice->TreeH();
162 Int_t ntracks = TH->GetEntries();
163 cout<<"ntracks "<<ntracks<<endl;
165 // Get pointers to Alice detectors and Digits containers
166 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
167 TClonesArray *Particles = gAlice->Particles();
170 // fill modules with sorted by module hits
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();
179 // get the Tree for clusters
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);
185 AliITSgeom *geom = ITS->GetITSgeom();
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;
198 if(SPDthickness < 200) {
199 ylim = SPDthickness/2 - 4;
201 ylim = SPDthickness/2 - 10;
204 for (Int_t idettype=0;idettype<3;idettype++) {
206 TClonesArray *ITSclusters = ITS->ClustersAddress(idettype);
207 //printf ("ITSclusters %p \n",ITSclusters);
209 if (idettype != 0) continue;
215 for (Int_t mod=0; mod<nent; mod++) {
216 AliITSmodule *itsModule = (AliITSmodule*)ITSmodules->At(mod);
217 geom->GetModuleId(mod,lay,lad,det);
219 Int_t nhits = itsModule->GetNhits();
220 //if(nhits) printf("module nhits %d %d\n",mod,nhits);
223 ITS->ResetClusters();
225 Int_t nclust = ITSclusters->GetEntries();
226 if (!nclust) continue;
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);
236 Int_t clustersizex = itsclu->NclX();
237 Int_t clustersizez = itsclu->NclZ();
238 Int_t xstart = itsclu->XStartf();
239 Int_t xstop = itsclu->XStopf();
240 Float_t fxstart = xstart*xpitch;
241 Float_t fxstop = (xstop+1)*xpitch;
242 Float_t zstart = itsclu->ZStart();
243 Float_t zstop = itsclu->ZStop();
244 Int_t zend = itsclu->Zend();
245 Int_t ntrover = itsclu->NTracks();
246 Float_t clusterx = 0;
247 Float_t clusterz = itsclu->Z();
250 itsclu->GetTracks(tr0,tr1,tr2);
252 for(Int_t ii=0;ii<clustersizex;ii++) {
253 clusterx += (xstart+0.5+ii)*xpitch;
255 clusterx /= clustersizex;
256 clusterz /= clustersizez;
258 Float_t clusterQ = itsclu->Q();
260 if(lay == 1) occup1 += clusterQ;
261 if(lay == 2) occup2 += clusterQ;
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;
273 Float_t dxprimlast = 10.e+6;
274 Float_t dzprimlast = 10.e+6;
276 Float_t xhit0 = 1e+6;
277 Float_t zhit0 = 1e+6;
278 for (Int_t hit=0;hit<nhits;hit++) {
280 // Find coordinate differences between the hit and cluster positions
281 // for the resolution determination.
283 itsHit = (AliITShit*)itsModule->GetHit(hit);
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();
291 Int_t hitstat = itsHit->GetTrackStatus();
292 Float_t zhit = 10000*itsHit->GetZL();
293 Float_t xhit = 10000*itsHit->GetXL();
294 Float_t yhit = 10000*itsHit->GetYL();
295 Int_t parent = itsHit->GetParticle()->GetFirstMother();
296 Int_t partcode = itsHit->GetParticle()->GetPdgCode();
300 if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery
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);
309 if(zhit > SPDlength/2) {
310 //cout<<"!!! z outside ="<<zhit<<endl;
311 zhit = SPDlength/2 - 10;
313 if(zhit < 0 && zhit < -SPDlength/2) {
314 //cout<<"!!! z outside ="<<zhit<<endl;
315 zhit = -SPDlength/2 + 10;
317 if(xhit > SPDwidth/2) {
318 //cout<<"!!! x outside ="<<xhit<<endl;
319 xhit = SPDwidth/2 - 10;
321 if(xhit < 0 && xhit < -SPDwidth/2) {
322 //cout<<"!!! x outside ="<<xhit<<endl;
323 xhit = -SPDwidth/2 + 10;
331 if(hitlayer == 1 && hitstat == 66 && yhit > ylim) {
332 //if(hitstat == 66 && hitprim ==1) {
338 if(hitlayer == 2 && hitstat == 66 && yhit < -ylim) {
343 if(hitstat != 68) continue; // Take only the hit if the last
344 // track point went out from
347 if(xhit0 > 1e+5 || zhit0 > 1e+5) continue;
349 Float_t xmed = (xhit + xhit0)/2;
350 Float_t zmed = (zhit + zhit0)/2;
352 Float_t xdif = xmed - clusterx;
353 Float_t zdif = zmed - clusterz;
357 // Consider the hits inside of cluster region only b.b.
359 if((xmed >= fxstart && xmed <= fxstop) && (zmed >= zstart && zmed <= zstop)) {
361 // Consider the hits only with the track number equaled to one
363 //if((track == tr0) || (track == tr1) || (track == tr2)) {
367 // part = (TParticle *)particles.UncheckedAt(track);
368 // Int_t partcode = part->GetPdgCode();
369 // Int_t primery = gAlice->GetPrimary(track);
372 // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
373 // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
375 Float_t pmod = itsHit->GetParticle()->P(); // total momentum at the
377 Float_t energy = itsHit->GetParticle()->Energy(); // energy at the
379 Float_t mass = itsHit->GetParticle()->GetMass(); // particle mass
381 Float_t pz = itsHit->GetParticle()->Pz(); // z momentum componetnt
383 Float_t px = itsHit->GetParticle()->Px(); // z momentum componetnt
385 Float_t py = itsHit->GetParticle()->Py(); // z momentum componetnt
387 Float_t phi = itsHit->GetParticle()->Phi(); // Phi angle at the
389 Float_t theta = itsHit->GetParticle()->Theta(); // Theta angle at the
391 //Float_t eta = itsHit->GetParticle()->Eta(); // Pseudo rapidity at the
394 if((energy-pz) > 0) {
395 y = 0.5*TMath::Log((energy+pz)/(energy-pz));
397 cout<<" Warning: energy < pz ="<<energy<<","<<pz<<endl;
400 Float_t eta = -TMath::Log(TMath::Tan(theta/2));
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);
411 if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
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
420 //if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery
423 if(hitprim > 0) noverprim = noverprim + 1;
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;
440 ntuple_st.pmod = pmod;
455 if(hitprim > 0) { // for primary particles
460 Ptot1->Fill(pmod/1000.);
472 Ptot2->Fill(pmod/1000.);
481 } // primery particles
483 } // end of cluster region
491 if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
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;
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;
526 // Write and Draw Histogramms and ntuples
530 TFile fhistos("SPD_his_dubna.root","RECREATE");
574 cout<<"!!! Histogramms and ntuples were written"<<endl;
576 TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
579 gPad->SetFillColor(33);
580 Xres1->SetFillColor(42);
583 gPad->SetFillColor(33);
584 Zres1->SetFillColor(46);
587 gPad->SetFillColor(33);
588 Xres2->SetFillColor(42);
591 gPad->SetFillColor(33);
592 Zres2->SetFillColor(46);
595 cout<<"END test for clusters and hits "<<endl;