3 void SPDclusterTest (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");
19 // Connect the Root Galice file containing Geometry, Kine and Hits
21 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
22 if (!file) file = new TFile("galice.root");
25 // Get AliRun object from file or create it if not on file
28 gAlice = (AliRun*)file->Get("gAlice");
29 if (gAlice) printf("AliRun object found on file\n");
30 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
38 for (int nev=0; nev<= evNumber2; nev++) {
39 Int_t nparticles = gAlice->GetEvent(nev);
40 cout << "nev " << nev <<endl;
41 cout << "nparticles " << nparticles <<endl;
42 if (nev < evNumber1) continue;
43 if (nparticles <= 0) return;
45 TTree *TH = gAlice->TreeH();
46 Int_t ntracks = TH->GetEntries();
47 cout<<"ntracks "<<ntracks<<endl;
51 AliITSRawClusterSPD *ITSclust;
53 // Get pointers to Alice detectors and Digits containers
54 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
55 TClonesArray *Particles = gAlice->Particles();
58 // fill modules with sorted by module hits
60 ITS->InitModules(-1,nmodules);
61 // ITS->FillModules(nev,-1,evNumber2,nmodules," "," ");
62 ITS->FillModules(nev,evNumber2,nmodules," "," ");
63 //get pointer to modules array
64 TObjArray *ITSmodules = ITS->GetModules();
67 // get the Tree for clusters
69 TTree *TC=ITS->TreeC();
70 Int_t nent=TC->GetEntries();
71 printf("Found %d entries in the tree (must be one per module per event!)\n",nent);
73 for (Int_t idettype=0;idettype<3;idettype++) {
75 TClonesArray *ITSclusters = ITS->ClustersAddress(idettype);
76 //printf ("ITSclusters %p \n",ITSclusters);
78 if (idettype != 0) continue;
81 // ------------ Cluster and point analysis histogramms ------------
83 TH1F *Nxpix1 = new TH1F("Nxpix1","Cluster size in x(r*phi) direction for layer 1",20,0.,20.);
84 TH1F *Nxpix2 = new TH1F("Nxpix2","Cluster size in x(r*phi) direction for layer 2",20,0.,20.);
85 TH1F *Nzpix1 = new TH1F("Nzpix1","Cluster size in z direction for layer 1",15,0.,15.);
86 TH1F *Nzpix2 = new TH1F("Nzpix2","Cluster size in z direction for layer 2",15,0.,15.);
87 TH1F *Xpix1 = new TH1F("Xpix1","Local x coordinate (mm) for layer 1",20,-2.,18.);
88 TH1F *Xpix2 = new TH1F("Xpix2","Local x coordinate (mm) for layer 2",20,-2.,18.);
89 TH1F *Zpix1 = new TH1F("Zpix1","Local z coordinate (mm) for layer 1",90,-2.,88.);
90 TH1F *Zpix2 = new TH1F("Zpix2","Lolac z coordinate (mm) for layer 2",90,-2.,88.);
92 TH1F *Xres1 = new TH1F("Xres1","Xrec and Xgen difference (micr) for layers 1",100,-200.,200.);
93 TH1F *Xres2 = new TH1F("Xres2","Xrec and Xgen difference (micr) for layers 2",100,-200.,200.);
94 TH1F *Zres1 = new TH1F("Zres1","Zrec and Zgen difference (micr) for layers 1",100,-800.,800.);
95 TH1F *Zres2 = new TH1F("Zres2","Zrec and Zgen difference (micr) for layers 2",100,-800.,800.);
98 // -------------- Create ntuples --------------------
100 // ntuple structures:
135 ntuple = new TTree("ntuple","Demo ntuple");
136 ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
137 ntuple->Branch("nx",&ntuple_st.nx,"nx/I");
138 ntuple->Branch("nz",&ntuple_st.nz,"nz/I");
139 ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
140 ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
141 ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
142 ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
143 ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");
145 ntuple1 = new TTree("ntuple1","Demo ntuple1");
146 ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
147 ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
148 ntuple1->Branch("det",&ntuple1_st.det,"det/I");
149 ntuple1->Branch("nx",&ntuple1_st.nx,"nx/I");
150 ntuple1->Branch("nz",&ntuple1_st.nz,"nz/I");
151 ntuple1->Branch("qcl",&ntuple1_st.qcl,"qcl/F");
152 ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I");
153 ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
154 ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
155 ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
156 ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
159 ntuple2 = new TTree("ntuple2","Demo ntuple2");
160 // ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I");
161 ntuple2->Branch("nx",&ntuple2_st.nx,"nx/I");
162 ntuple2->Branch("nz",&ntuple2_st.nz,"nz/I");
164 // ------------------------------------------------------------------------
168 for (Int_t mod=0; mod<nent; mod++) {
169 AliITSmodule *itsModule = (AliITSmodule*)ITSmodules->At(mod);
171 Int_t nhits = itsModule->GetNhits();
172 if(nhits) printf("module nhits %d %d\n",mod,nhits);
175 ITS->ResetClusters();
177 Int_t nclust = ITSclusters->GetEntries();
178 if (!nclust) continue;
182 for (Int_t clu=0;clu<nclust;clu++) {
183 itsclu = (AliITSRawClusterSPD*)ITSclusters->UncheckedAt(clu);
188 Int_t clustersizex = itsclu->NclX();
189 Int_t clustersizez = itsclu->NclZ();
190 // Int_t xstart = itsclu->XStart();
191 // Int_t xstop = itsclu->XStop();
192 Int_t xstart = itsclu->XStartf();
193 Int_t xstop = itsclu->XStopf();
194 Float_t fxstart = xstart*50;
195 Float_t fxstop = (xstop+1)*50;
196 Float_t zstart = itsclu->ZStart();
197 Float_t zstop = itsclu->ZStop();
198 Int_t zend = itsclu->Zend();
199 Int_t ntrover = itsclu->NTracks();
200 Float_t clusterx = itsclu->X();
201 Float_t clusterz = itsclu->Z();
202 Float_t clusterQ = itsclu->Q();
204 ntuple2_st.nx = clustersizex;
205 ntuple2_st.nz = clustersizez;
210 Float_t dxprimlast = 10.e+6;
211 Float_t dzprimlast = 10.e+6;
213 if(clustersizex>2&&clustersizez>1) {
214 // if(module > 217 && module < 226) {
215 if (nclust) printf("Found %d clust for module %d in det type %d \n",nclust,mod,idettype);
216 cout<<"mod,nclust,clu,Nxpix,Nzpix ="<<mod<<","<<nclust<<","<<clu<<","<<clustersizex<<","<<clustersizez<<endl;
217 // cout<<"clusx,clusz ="<<clusterx<<","<<clusterz<<endl;
218 cout<<"XStartf,XStopf,Zend,Zstart,Zstop,Q ="<<xstart<<","<<xstop<<","<<zend<<","<<zstart<<","<<zstop<<","<<clusterQ<<endl;
222 Float_t SPDlength = 83600;
223 Float_t SPDwidth = 12800;
224 Float_t xhit0 = 1e+5;
225 Float_t zhit0 = 1e+5;
228 for (Int_t hit=0;hit<nhits;hit++) {
230 // Find coordinate differences between the hit and cluster positions
231 // for the resolution determination.
233 itsHit = (AliITShit*)itsModule->GetHit(hit);
235 Int_t hitlayer = itsHit->GetLayer();
236 Int_t hitladder= itsHit->GetLadder();
237 Int_t hitdet= itsHit->GetDetector();
239 Int_t clusterlayer = hitlayer;
240 Int_t clusterladder= hitladder;
241 Int_t clusterdetector = hitdet;
243 Int_t track = itsHit->fTrack;
245 Int_t hitstat = itsHit->GetTrackStatus();
248 Float_t zhit = 10000*itsHit->GetZL();
249 Float_t xhit = 10000*itsHit->GetXL();
251 if(abs(zhit) > SPDlength/2) {
252 if(hitstat == 66) zhit0 = 1e+5;
256 if(abs(xhit) > SPDwidth/2) {
257 if(hitstat == 66) xhit0 = 1e+5;
263 Float_t yhit = 10000*itsHit->GetYL();
265 if(hitlayer == 1 && hitstat == 66 && yhit > 71) {
269 if(hitlayer == 2 && hitstat == 66 && yhit < -71) {
274 if(hitstat != 68) continue; // Take only the hit if the last
275 // track point went out from
278 if(xhit0 > 9e+4 || zhit0 > 9e+4) continue;
281 Float_t xmed = (xhit + xhit0)/2;
282 Float_t zmed = (zhit + zhit0)/2;
284 Float_t xdif = xmed - clusterx;
285 Float_t zdif = zmed - clusterz;
287 // cout<<"clu,hit,xmed,fxstart,fxstop,zmed,zstart,zstop ="<<clu<<","<<hit<<","<<xmed<<","<<fxstart<<","<<fxstop<<","<<zmed<<","<<zstart<<","<<zstop<<endl;
289 // Consider the hits inside of cluster region only
291 if((xmed >= fxstart && xmed <= fxstop) && (zmed >= zstart && zmed <= zstop)) {
295 // part = (TParticle *)particles.UncheckedAt(track);
296 // Int_t partcode = part->GetPdgCode();
297 // Int_t primery = gAlice->GetPrimary(track);
299 Int_t parent = itsHit->GetParticle()->GetFirstMother();
300 Int_t partcode = itsHit->GetParticle()->GetPdgCode();
302 // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
303 // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
306 Float_t px = itsHit->GetPXL();
307 Float_t py = itsHit->GetPYL();
308 Float_t pz = itsHit->GetPZL();
309 Float_t pmod = 1000*sqrt(px*px+py*py+pz*pz);
311 // cout<<"track,partcode,pmod,parent ="<<track<<","<<partcode<<","<<pmod<<","<<parent<<endl;
315 if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
318 if(dray == 0) noverlaps = noverlaps + 1; // overlapps for all hits but
319 // not for delta ray which
320 // also went out from the
321 // detector and returned
324 if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery
327 if(hitprim > 0) noverprim = noverprim + 1;
336 ntuple_st.lay = hitlayer;
337 ntuple_st.nx = clustersizex;
338 ntuple_st.nz = clustersizez;
339 ntuple_st.hitprim = hitprim;
340 ntuple_st.partcode = partcode;
343 ntuple_st.pmod = pmod;
349 if(hitprim > 0) { // for primary particles
351 // cout<<"!!!!!! lay,hitprim,xdif,zdif ="<<hitlayer<<","<<hitprim<<","<<xdif<<","<<zdif<<endl;
356 // cout<<"!!!!!! lay,hitprim,xdif,zdif ="<<hitlayer<<","<<hitprim<<","<<xdif<<","<<zdif<<endl;
360 } // primery particles
362 } // end of cluster region
370 // ntuple1->Fill(clusterlayer,clustersizex,clustersizez,noverlaps,\
373 if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
376 ntuple1_st.lay = clusterlayer;
377 ntuple1_st.lad = clusterladder;
378 ntuple1_st.det = clusterdetector;
379 ntuple1_st.nx = clustersizex;
380 ntuple1_st.nz = clustersizez;
381 ntuple1_st.qcl = clusterQ;
382 ntuple1_st.ntrover = ntrover;
383 ntuple1_st.noverlaps = noverlaps;
384 ntuple1_st.noverprim = noverprim;
385 ntuple1_st.dx = dxprimlast;
386 ntuple1_st.dz = dzprimlast;
398 // Write and Draw Histogramms and ntuples
402 TFile fhistos("SPD_his.root","RECREATE");
424 cout<<"!!! Histogramms and ntuples were written"<<endl;
426 TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
429 gPad->SetFillColor(33);
430 Xres1->SetFillColor(42);
433 gPad->SetFillColor(33);
434 Zres1->SetFillColor(46);
437 gPad->SetFillColor(33);
438 Xres2->SetFillColor(42);
441 gPad->SetFillColor(33);
442 Zres2->SetFillColor(46);
445 cout<<"END test for clusters and hits "<<endl;