EMCAL/DCAL Trigger Mapping for Run 2
[u/mrichter/AliRoot.git] / ITS / ShowSPDRecPoints.C
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
51 void 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){                        
206         delete AliRunLoader::Instance();   
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");
216       man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
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 }