]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/ShowSPDRecPoints.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix bug in pileup tagging (F. Prino)
[u/mrichter/AliRoot.git] / ITS / ShowSPDRecPoints.C
CommitLineData
d55d87f0 1//////////////////////////////////////////////////////////////////
2// Macro to check clusters in the 2 SPD layers //
3// Provides: //
4// 2 canvases with //
5// - cluster loc and glob coordinates (each layer) //
6// 1 canvas with //
7// - cluster glob coordinates 2D and 3D //
8// 1 canvas with //
9// - correlations of #clusters for sectors //
10// - correlations of #clusters for half-sectors //
11// //
12// Maria.Nicassio@ba.infn.it //
13// Domenico.Elia@ba.infn.it //
14//////////////////////////////////////////////////////////////////
15
16#if !defined(__CINT__) || defined(__MAKECINT__)
17
18#include <Riostream.h>
19#include <TFile.h>
20#include <TTree.h>
21#include <TBranch.h>
22#include <TCanvas.h>
23#include <TClonesArray.h>
24#include <TH1.h>
25#include <TH2.h>
26#include <TH3.h>
27#include <TStyle.h>
28#include <TProfile.h>
29#include <TClassTable.h>
30#include <TInterpreter.h>
31#include <TGraph2D.h>
32#include <TGraph.h>
33#include <TGeoManager.h>
34
35#include "AliCDBManager.h"
36#include "AliRunLoader.h"
37#include "AliESD.h"
38#include "AliRun.h"
39#include "AliGeomManager.h"
40#include "AliITS.h"
41#include "AliITSgeom.h"
42#include "AliITSLoader.h"
43#include <AliITSRecPoint.h>
44#include <AliHeader.h>
45#include <AliITSDetTypeRec.h>
46
47#endif
48
49/* $Id$ */
50
51void ShowSPDRecPoints(Int_t RunStart, Int_t RunStop){
52
53 // Set data directory
54 Char_t str[256];
55 Char_t* dir = "/data/alipix/PhysicsEvents/pp/ppProd/pdc07/oldgeom";
56
57 // Variables for histo booking and filling
58 Int_t modmin=0;
59 Int_t modmax=240;
60 Int_t totmod=modmax-modmin;
61
62 Float_t xlim[2]={4.5,8.};
63 Float_t zlim[2]={15.,15.};
64
65 Int_t nClusters[2];
66
67 Float_t clustersXCoord[2][200];
68 Float_t clustersYCoord[2][200];
69 Float_t clustersZCoord[2][200];
70
71 Float_t cluGlo[3]={0.,0.,0.};
72
73 Int_t nClustersPerLayerPerSector[2][10];
74 Int_t nClustersPerLayerPerHalfSector[2][20];
75
76 Int_t iSector=0;
77 Int_t iHalfSector=0;
78
79 // Booking of histograms
80 gStyle->SetPalette(1,0);
81 TH1F* hlayer= new TH1F("hlayer","",6,0.,6.);
82 TH1F** hmod = new TH1F*[2];
83 TH1F** hxl = new TH1F*[2];
84 TH1F** hzl = new TH1F*[2];
85 TH1F** hxg = new TH1F*[2];
86 TH1F** hyg = new TH1F*[2];
87 TH1F** hzg = new TH1F*[2];
88 TH1F** hr = new TH1F*[2];
89 TH1F** hphi = new TH1F*[2];
90 TH1F** hMultSPDcl = new TH1F*[2];
91 TH1F** hType = new TH1F*[2]; // cluster type ?
92
93 TH2F** hr_phi = new TH2F*[2];
94 TH2F** hx_y = new TH2F*[2];
95 TH3F** hx_y_z = new TH3F*[2];
96
97 TH2F** hnCl2_nCl1_Sec = new TH2F*[10];
98 TProfile** hnCl2vsnCl1_Sec = new TProfile*[10];
99 TH2F** hnCl2_nCl1_HSec = new TH2F*[20];
100// TProfile** hnCl2vsnCl1_HSec = new TProfile*[20];
101 TH2F** hNy_Nz = new TH2F*[2]; // y and z length
102 TH2F** hPhi_Z = new TH2F*[2];
103
104 TH2F *hMultSPDcl2_MultSPDcl1 =
105 new TH2F("nCl2_nCl1","# SPD clusters on layer 2 vs # on layer 1",200,0.,200.,200,0.,200.);
106 hMultSPDcl2_MultSPDcl1->GetXaxis()->SetTitle("# clusters (inner layer)");
107 hMultSPDcl2_MultSPDcl1->GetYaxis()->SetTitle("# clusters (outer layer)");
108 Char_t name[10];
109 Char_t title[20];
110 for (Int_t iLay=0;iLay<2;iLay++) {
111 sprintf(name,"hmod%d",iLay+1);
112 hmod[iLay]=new TH1F(name,"SPD clusters - Module number",totmod,modmin,modmax);
113 hmod[iLay]->GetXaxis()->SetTitle("Module number");
114 hmod[iLay]->GetYaxis()->SetTitle("Entries");
115 sprintf(name,"hxloc%d",iLay+1);
116 hxl[iLay]=new TH1F(name,"SPD clusters - Local x coordinate",100,-4.,4.);
117 hxl[iLay]->GetXaxis()->SetTitle("Local x [cm]");
118 hxl[iLay]->GetYaxis()->SetTitle("Entries");
119 sprintf(name,"hzloc%d",iLay+1);
120 hzl[iLay]=new TH1F(name,"SPD clusters - Local z coordinate",100,-4.,4.);
121 hzl[iLay]->GetXaxis()->SetTitle("Local z [cm]");
122 hzl[iLay]->GetYaxis()->SetTitle("Entries");
123 sprintf(name,"hxgl%d",iLay+1);
124 hxg[iLay]=new TH1F(name,"SPD clusters - Global x coordinate",100,-xlim[iLay],xlim[iLay]);
125 hxg[iLay]->GetXaxis()->SetTitle("Global x [cm]");
126 hxg[iLay]->GetYaxis()->SetTitle("Entries");
127 sprintf(name,"hygl%d",iLay+1);
128 hyg[iLay]=new TH1F(name,"SPD clusters - Global y coordinate",100,-xlim[iLay],xlim[iLay]);
129 hyg[iLay]->GetXaxis()->SetTitle("Global y [cm]");
130 hyg[iLay]->GetYaxis()->SetTitle("Entries");
131 sprintf(name,"hzgl%d",iLay+1);
132 hzg[iLay]=new TH1F(name,"SPD clusters - Global z coordinate",150,-zlim[iLay],zlim[iLay]);
133 hzg[iLay]->GetXaxis()->SetTitle("Global z [cm]");
134 hzg[iLay]->GetYaxis()->SetTitle("Entries");
135 sprintf(name,"hr%d",iLay+1);
136 hr[iLay]=new TH1F(name,"SPD clusters - r",100,0.,50.);
137 hr[iLay]->GetXaxis()->SetTitle("r [cm]");
138 hr[iLay]->GetYaxis()->SetTitle("Entries");
139 sprintf(name,"hphi%d",iLay+1);
140 hphi[iLay]=new TH1F(name,"SPD clusters - #phi",100,0.,2*TMath::Pi());
141 hphi[iLay]->GetXaxis()->SetTitle("#phi [rad]");
142 hphi[iLay]->GetYaxis()->SetTitle("Entries");
143 sprintf(name,"hType%d",iLay+1);
144 hType[iLay]=new TH1F(name,"SPD clusters - Type",100,0.,100.);
145 hType[iLay]->GetXaxis()->SetTitle("Cluster type");
146 hType[iLay]->GetYaxis()->SetTitle("Entries");
147
148 sprintf(name,"hMultSPDcl%d",iLay+1);
149 hMultSPDcl[iLay]=new TH1F(name,"Cluster multiplicity",200,0.,200.);
150 hMultSPDcl[iLay]->GetXaxis()->SetTitle("Cluster multiplicity");
151 hMultSPDcl[iLay]->GetYaxis()->SetTitle("Entries");
152
153 sprintf(name,"hNy_Nz%d",iLay+1);
154 hNy_Nz[iLay]=new TH2F(name,"SPD clusters - Length",100,0.,100.,100,0.,100.);
155 hNy_Nz[iLay]->GetXaxis()->SetTitle("z length");
156 hNy_Nz[iLay]->GetYaxis()->SetTitle("y length");
157 sprintf(name,"hPhi_Z%d",iLay+1);
158 hPhi_Z[iLay]=new TH2F(name,"SPD clusters - #phi vs z",150,-zlim[iLay],zlim[iLay],100,0.,2*TMath::Pi());
159 hPhi_Z[iLay]->GetXaxis()->SetTitle("Global z [cm]");
160 hPhi_Z[iLay]->GetYaxis()->SetTitle("#phi [rad]");
161 sprintf(name,"hr_phi%d",iLay+1);
162 hr_phi[iLay]=new TH2F(name,"SPD clusters - #phi_r",500,0.,50.,100,0.,2*TMath::Pi());
163 hr_phi[iLay]->GetXaxis()->SetTitle("r [cm]");
164 hr_phi[iLay]->GetYaxis()->SetTitle("#phi [rad]");
165 sprintf(name,"hx_y%d",iLay+1);
166 hx_y[iLay]=new TH2F(name,"SPD clusters - y_x",200,-10.,10.,200,-10.,10.);
167 hx_y[iLay]->GetXaxis()->SetTitle("x [cm]");
168 hx_y[iLay]->GetYaxis()->SetTitle("y [cm]");
169 sprintf(name,"hx_y_z%d",iLay+1);
170 hx_y_z[iLay]=new TH3F(name,"SPD clusters - y_x",200,-10.,10.,200,-10.,10.,150,-15.,15.);
171 hx_y_z[iLay]->GetXaxis()->SetTitle("z [cm]");
172 hx_y_z[iLay]->GetYaxis()->SetTitle("x [cm]");
173 hx_y_z[iLay]->GetZaxis()->SetTitle("y [cm]");
174 }
175 for (Int_t iSec=0; iSec<10; iSec++) {
176 sprintf(name,"hnCl2_nCl1_Sector%d",iSec);
177 sprintf(title,"Sector %d",iSec+1);
178 hnCl2_nCl1_Sec[iSec]=new TH2F(name,title,200,0.,200.,200,0.,200.);
179 hnCl2_nCl1_Sec[iSec]->GetXaxis()->SetTitle("# clusters (inner layer)");
180 hnCl2_nCl1_Sec[iSec]->GetYaxis()->SetTitle("# clusters (outer layer)");
181 sprintf(name,"hnCl2vsnCl1_Sector%d",iSec);
182 sprintf(title,"Sector %d",iSec+1);
183 hnCl2vsnCl1_Sec[iSec]=new TProfile(name,title,200,0.,200.,0.,200.);
184 hnCl2vsnCl1_Sec[iSec]->GetXaxis()->SetTitle("# clusters (inner layer)");
185 hnCl2vsnCl1_Sec[iSec]->GetYaxis()->SetTitle("# clusters (outer layer)");
186 }
187 for (Int_t iHalfSec=0; iHalfSec<20; iHalfSec++) {
188 sprintf(name,"hnCl2_nCl1_HalfSector%d",iHalfSec);
189 sprintf(title,"Half-Sector %d",iHalfSec+1);
190 hnCl2_nCl1_HSec[iHalfSec]=new TH2F(name,title,200,0.,200.,200,0.,200.);
191 hnCl2_nCl1_HSec[iHalfSec]->GetXaxis()->SetTitle("# clusters (inner layer)");
192 hnCl2_nCl1_HSec[iHalfSec]->GetYaxis()->SetTitle("# clusters (outer layer)");
193 }
194
195 // Loop over runs
196 for (Int_t run=RunStart; run<RunStop+1; run++) {
197
198 // setup galice and runloader
199// cout << "File nr --> " << run << endl;
200
201 if (gClassTable->GetID("AliRun") < 0) {
202 gInterpreter->ExecuteMacro("loadlibs.C");
203 }
204 else {
205 if (gAlice){
33c3c91a 206 delete AliRunLoader::Instance();
d55d87f0 207 delete gAlice;
208 gAlice=0;
209 }
210 }
211
212 // Set OfflineConditionsDataBase if needed
213 AliCDBManager* man = AliCDBManager::Instance();
214 if (!man->IsDefaultStorageSet()) {
215 printf("Setting a local default storage and run number 0\n");
162637e4 216 man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
d55d87f0 217 man->SetRun(0);
218 }
219 else {
220 printf("Using deafult storage \n");
221 }
222
223 // retrives geometry
224 if (!gGeoManager) {
225 sprintf(str,"%s/ppMinBias%04d/geometry.root",dir,run);
226 AliGeomManager::LoadGeometry(str);
227 }
228
229 sprintf(str,"%s/ppMinBias%04d/galice.root",dir,run);
230 AliRunLoader* rl = AliRunLoader::Open(str);
231
232 if (rl == 0x0){
233 cerr<<"Can not open session RL=NULL"<< endl;
234 return;
235 }
236 Int_t retval = rl->LoadgAlice();
237 if (retval){
238 cerr<<"LoadgAlice returned error"<<endl;
239 return;
240 }
241 gAlice=rl->GetAliRun();
242
243 retval = rl->LoadHeader();
244 if (retval){
245 cerr<<"LoadHeader returned error"<<endl;
246 return;
247 }
248
249 AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
250 if (!ITSloader){
251 cerr<<"ITS loader not found"<<endl;
252 return;
253 }
254 ITSloader->LoadRecPoints("read");
255
256 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
257 ITS->SetTreeAddress();
258 AliITSgeom *geom = ITS->GetITSgeom();
259 if (!geom) {
260 cout << " Can't get the ITS geometry!" << endl;
261 return ;
262 }
263
264 AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
265
266 detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
267 detTypeRec->SetDefaults();
268
269 Int_t nEvents=rl->GetNumberOfEvents();
270 printf("Total Number of events = %d\n",nEvents);
271
272 for (Int_t iev=0;iev<nEvents;iev++){
273 rl->GetEvent(iev);
274
275 TTree *TR = ITSloader->TreeR();
276 TClonesArray* ITSClusters;
277 ITSClusters = detTypeRec->RecPoints();
278 TBranch *branch = 0;
279 if (TR && ITSClusters) {
280 branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
281 if (branch) branch->SetAddress(&ITSClusters);
282 }
283
284// Int_t nparticles = rl->GetHeader()->GetNtrack();
285// cout<<"Event # "<<iev<<" #Particles="<<nparticles<<endl;
286
287 // Reset cluster counters
288 for (Int_t iLay=0; iLay<2; iLay++) {
289 nClusters[iLay]=0;
290 for (Int_t iSec=0; iSec<10; iSec++) {
291 nClustersPerLayerPerSector[iLay][iSec]=0;
292 }
293 for (Int_t iHSec=0; iHSec<20; iHSec++) {
294 nClustersPerLayerPerHalfSector[iLay][iHSec]=0;
295 }
296 }
297 for (Int_t subd=0;subd<1;subd++) {
298
299 Int_t first = geom->GetStartDet(subd);
300 Int_t last = geom->GetLastDet(subd);
301 for (Int_t mod=first; mod<=last; mod++) {
302 detTypeRec->ResetRecPoints();
303 branch->GetEvent(mod);
304 Int_t nrecp = ITSClusters->GetEntries();
305 while (nrecp--) {
306
307 AliITSRecPoint *recp = (AliITSRecPoint*)ITSClusters->At(nrecp);
308
309 Int_t lay=recp->GetLayer();
310// cout<<"lay"<<lay<<endl;
311 recp->GetGlobalXYZ(cluGlo);
312
313 Float_t rad=TMath::Sqrt(cluGlo[0]*cluGlo[0]+cluGlo[1]*cluGlo[1]);
314 Float_t phi= TMath::Pi() + TMath::ATan2(-cluGlo[1],-cluGlo[0]);
315 hlayer->Fill(lay);
316 hmod[lay]->Fill(mod);
317 hzl[lay]->Fill(recp->GetDetLocalZ());
318 hxl[lay]->Fill(recp->GetDetLocalX());
319 hzg[lay]->Fill(cluGlo[2]);
320 hyg[lay]->Fill(cluGlo[1]);
321 hxg[lay]->Fill(cluGlo[0]);
322 hr[lay]->Fill(rad);
323 hphi[lay]->Fill(phi);
324 hPhi_Z[lay]->Fill(cluGlo[2],phi);
325 hr_phi[lay]->Fill(rad,phi);
326 hx_y[lay]->Fill(cluGlo[0],cluGlo[1]);
327 hx_y_z[lay]->Fill(cluGlo[2],cluGlo[0],cluGlo[1]);
328 clustersXCoord[lay][nClusters[lay]]=cluGlo[0];
329 clustersYCoord[lay][nClusters[lay]]=cluGlo[1];
330 clustersZCoord[lay][nClusters[lay]]=cluGlo[2];
331 nClusters[lay]++;
332
333 hType[lay]->Fill(recp->GetType());
334// cout<<"Clusters type"<<recp->GetType()<<endl;
335 hNy_Nz[lay]->Fill(recp->GetNz(),recp->GetNy());
336
337 // Set Sector number and increase the counter
338 if (lay==0) {
339 for (Int_t nRange=0; nRange<10; nRange++) {
340 if ((mod>=nRange*8) && (mod<=(nRange*8+7))) {
341 iSector = nRange;
342 nClustersPerLayerPerSector[lay][iSector]++;
343 break;
344 }
345 }
346 }
347 if (lay==1) {
348 for (Int_t nRange=0; nRange<10; nRange++) {
349 if ((mod>=80+nRange*16) && (mod<=(80+nRange*16+15))) {
350 iSector = nRange;
351 nClustersPerLayerPerSector[lay][iSector]++;
352 break;
353 }
354 }
355 }
356 // Set HalfSector number and increase the counter
357 if (lay==0) {
358 for (Int_t nRange=0; nRange<20; nRange++) {
359 if ((mod>=nRange*4) && (mod<=(nRange*4+3))) {
360 iHalfSector = nRange;
361 nClustersPerLayerPerHalfSector[lay][iHalfSector]++;
362 break;
363 }
364 }
365 } else {
366 for (Int_t nRange=0; nRange<20; nRange++) {
367 if ((mod>=80+nRange*4) && (mod<=(80+nRange*4+3))) {
368 iHalfSector = nRange;
369 nClustersPerLayerPerHalfSector[lay][iHalfSector]++;
370 break;
371 }
372 }
373 }
374 }
375 }
376 }
377
378 for (Int_t iLay=0; iLay<2; iLay++) hMultSPDcl[iLay]->Fill(nClusters[iLay]);
379 for (Int_t iSec=0; iSec<10; iSec++) {
380 hnCl2_nCl1_Sec[iSec]->Fill(nClustersPerLayerPerSector[0][iSec],nClustersPerLayerPerSector[1][iSec]);
381 hnCl2vsnCl1_Sec[iSec]->Fill(nClustersPerLayerPerSector[0][iSec],nClustersPerLayerPerSector[1][iSec]);
382 }
383 for (Int_t iHSec=0; iHSec<20; iHSec++) {
384 hnCl2_nCl1_HSec[iHSec]->Fill(nClustersPerLayerPerHalfSector[0][iHSec],nClustersPerLayerPerHalfSector[1][iHSec]);
385 }
386 hMultSPDcl2_MultSPDcl1->Fill(nClusters[0],nClusters[1]);
387
388 } //end loop over events
389 rl->UnloadAll();
390 delete rl;
391
392 } //end loop over runs
393
394 // Draw and Write histos
395
396 TFile* fout = new TFile("out_ShowSPDRecPoints.root","RECREATE");
397
398 hlayer->Write();
399// cev0->Write();
400
401 TCanvas **c=new TCanvas*[2];
402 Char_t ctit[30];
403 for(Int_t iLay=0;iLay<2;iLay++){
404 hNy_Nz[iLay]->Write();
405 sprintf(name,"can%d",iLay+1);
406 sprintf(ctit,"Layer %d",iLay+1);
407 c[iLay]=new TCanvas(name,ctit,1200,900);
408 c[iLay]->Divide(3,3);
409 c[iLay]->cd(1);
410 hmod[iLay]->Draw();
411 hmod[iLay]->Write();
412 c[iLay]->cd(2);
413 hxl[iLay]->Draw();
414 hxl[iLay]->Write();
415 c[iLay]->cd(3);
416 hzl[iLay]->Draw();
417 hzl[iLay]->Write();
418 c[iLay]->cd(4);
419 hxg[iLay]->Draw();
420 hxg[iLay]->Write();
421 c[iLay]->cd(5);
422 hyg[iLay]->Draw();
423 hyg[iLay]->Write();
424 c[iLay]->cd(6);
425 hzg[iLay]->Draw();
426 hzg[iLay]->Write();
427 c[iLay]->cd(7);
428 hr[iLay]->Draw();
429 hr[iLay]->Write();
430 c[iLay]->cd(8);
431 hphi[iLay]->Draw();
432 hphi[iLay]->Write();
433 c[iLay]->cd(9);
434 hPhi_Z[iLay]->Draw("colz");
435 hPhi_Z[iLay]->Write();
436 hr_phi[iLay]->Write();
437 hx_y[iLay]->Write();
438 hx_y_z[iLay]->Write();
439 }
440
441 TCanvas *cCoord=new TCanvas("cCoord","Cluster coordinates",1200,900);
442 cCoord->Divide(2,2);
443
444 for (Int_t iLay=0;iLay<2;iLay++) {
445 cCoord->cd(1);
446 hx_y[iLay]->SetMarkerStyle(22);
447 hx_y[iLay]->SetMarkerSize(0.3);
448 hx_y[iLay]->SetMarkerColor(iLay+1);
449 if (iLay==0) hx_y[iLay]->Draw("p");
450 else hx_y[iLay]->Draw("p,same");
451 cCoord->cd(2);
452 hx_y_z[iLay]->SetMarkerStyle(23);
453 hx_y_z[iLay]->SetMarkerSize(0.3);
454 hx_y_z[iLay]->SetMarkerColor(iLay+1);
455 hx_y_z[iLay]->Draw("p,same");
456 cCoord->cd(3);
457 hr_phi[iLay]->SetMarkerColor(iLay+1);
458 hr_phi[iLay]->Draw("p,same");
459 }
460 TCanvas *cCorrelations_Sectors=new TCanvas("cCorrelations_S","SPD cluster correlations (Sectors)",1200,900);
461 cCorrelations_Sectors->Divide(3,5);
462 cCorrelations_Sectors->cd(1);
463 hMultSPDcl2_MultSPDcl1->Draw();
464 hMultSPDcl2_MultSPDcl1->Write();
465 cCorrelations_Sectors->cd(2);
466 hMultSPDcl[0]->Draw();
467 hMultSPDcl[0]->Write();
468 cCorrelations_Sectors->cd(3);
469 hMultSPDcl[1]->Draw();
470 hMultSPDcl[1]->Write();
471 for (Int_t iS=0; iS<10; ++iS) {
472 cCorrelations_Sectors->cd(iS+4);
473 hnCl2_nCl1_Sec[iS]->Draw();
474 hnCl2_nCl1_Sec[iS]->Write();
475 hnCl2vsnCl1_Sec[iS]->Write();
476 }
477
478 TCanvas *cCorrelations_HSectors=new TCanvas("cCorrelations_HS","SPD cluster correlations (Half-Sectors)",1200,900);
479 cCorrelations_HSectors->Divide(5,4);
480
481 for (Int_t iHS=0; iHS<20; ++iHS) {
482 cCorrelations_HSectors->cd(iHS+1);
483 hnCl2_nCl1_HSec[iHS]->Write();
484 hnCl2_nCl1_HSec[iHS]->Draw();
485 }
486
487 fout->Close();
488
489}