]> git.uio.no Git - u/mrichter/AliRoot.git/blob - ITS/AliITSplotSPDClusters.C
Macros to use Dubna SPD simulations instead of the default Bari/Salerno,
[u/mrichter/AliRoot.git] / ITS / AliITSplotSPDClusters.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2
3 // Standard includes
4 #include <iostream.h>
5 // Root includes
6 #include <TFile.h>
7 #include <TObjArray.h>
8 #include <TClonesArray.h>
9 #include <TTree.h>
10 #include <TBranch.h>
11 #include <TH1D.h>
12 // AliRoot includes
13 #include "STEER/AliRun.h"
14 #include "ITS/AliITS.h"
15 #include "ITS/AliITSRawCluster.h"
16
17 #endif
18
19 void AliITSplotSPDClusters(const char* filename="galice_80.root"){
20     //------------------------------------------------------------------
21     // This macro will read the TreeC produced during reconstruction and
22     // plot the number and type of clusters found in the SPD.
23     // root [0] .L AliITSplotSPDClusters.C  //this loads the macro in memory
24     // root [1] AliITSplotSPDClusters();    //by default process first event
25     // or
26     // root [1] AliITSplotSPDClusters("galice2.root"); // process all events
27     //                                                    from the file
28     //                                                    galice2.root.
29     // or
30     // root [0] .x AliITSplotSPDClusters.C("galice2.root") // will all of the
31     //                                                        events from the
32     //                                                        galice2.root
33     //------------------------------------------------------------------
34
35     // delete the current gAlice object, the one from input file will be used
36     if(gAlice){
37         delete gAlice;
38         gAlice = 0;
39     }else{ // Dynamically link some shared libs
40         if(gClassTable->GetID("AliRun") < 0) {
41             gROOT->LoadMacro("loadlibs.C");
42             loadlibs();
43         } // end if
44     } // end if gAlice
45
46     // Connect the Root Galice file containing Geometry, Kine and Hits
47     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
48     if(!file) file = new TFile(filename);
49     // Get AliRun object from file or create it if not on file
50     if(!gAlice) {
51         gAlice = (AliRun*)file->Get("gAlice");
52         if(gAlice) cout << "AliRun object found on file" << endl;
53         if(!gAlice) gAlice = new AliRun("gAlice","SPD Clulster analysis");
54     } // end if !gAlice
55     Int_t nevents = gAlice->GetEventsPerRun();
56
57     // Get pointers to the ITS Alice detectors.
58     AliITS *ITS = (AliITS*) gAlice->GetDetector("ITS");
59     if(ITS==0){
60         cout << "ITS not found. Exiting." << endl;
61         return;
62     } // end if ITS==0
63     AliITSresponseSPDdubna *resp0 = ITS->DetType(0)->GetResponseModel();
64     Float_t diffCoeff= resp0->DistanceOverVoltage();//Get Value of Diffusion Coefficient parameter d/v.
65
66     TH1D *hclx1 = new TH1D("hclx1","Number of SPD Clusters in x, layer 1",
67                            15,0.5,15.5);
68     TH1D *hclz1 = new TH1D("hclz1","Number of SPD Clusters in z, layer 1",
69                            5,0.5,5.5);
70     TH1D *hclx2 = new TH1D("hclx2","Number of SPD Clusters in x, layer 2",
71                            15,0.5,15.5);
72     TH1D *hclz2 = new TH1D("hclz2","Number of SPD Clusters in z, layer 2",
73                            5,0.5,5.5);
74     // Create Arrays with clusters from:  Data, Ba/Sa Model, old version of
75     // Dubna
76     Float_t dataX1[9] = {0.493, 0.493, 0.0273, 0.00617, 0.00112, 0.000257,
77                          0.000257,  0.000257, 0.000257};
78                         //clusters from the data in the x-dimension with
79                         //standard paramaters
80     Float_t baSaX1[11] = {0.942, 0.0180, 0.0112, 0.00389, 0.00203, 0.000257,
81                           0.001, 0.000257, 0.000257, 0.000257, 0.001};
82                         //same for Ba/Sa model
83     Float_t dubnaX1[7] =  {0.889, 0.0789, 0.0151, 0.00389, 0.001, 0.001,
84                            0.001};
85                         //same for old version of Dubna model
86     Float_t dataX2[9] = {0.482, 0.482, 0.0325, 0.00791, 0.00157, 0.000275,
87                          0.000275, 0.000275, 0.000275};
88                         //clusters from the data in the x-dimension with
89                         //optimized paramaters
90     Float_t baSaX2[11] = {0.946, 0.0175, 0.0245, 0.00482, 0.001, 0.001,
91                           0.000275,  0.001, 0.000275, 0.000275, 0.001};
92                         //same for Ba/Sa model
93     Float_t dubnaX2[8] = {0.638, 0.325, 0.0275, 0.01, 0.00431, 0.000275,
94                           0.001, 0.001};
95                         //same for old version of Dubna model
96     Float_t dataZ1[4] = {0.889, 0.0624, 0.000257, 0.000257};
97                         //clusters from the data in the z-dimension with
98                         //standard paramaters
99     Float_t baSaZ1[3] = {1., 0.0160, 0.000257}; //same for Ba/Sa model
100     Float_t dubnaZ1[3] = {0.889, 0.0180, 0.000257};
101                         //same for old version of Dubna model
102     Float_t dataZ2[4] = {0.889, 0.0624, 0.000257, 0.000257};
103                         //clusters from the data in the z-dimension with
104                         //optimized paramaters
105     Float_t baSaZ2[4] = {1., 0.0160, 0.00215, 0.000257}; //same for Ba/Sa model
106     Float_t dubnaZ2[3] = {0.889, 0.0412, 0.000257};
107                         //same for old version of Dubna model
108     
109     TH1F *hdataX1  = new TH1F("hdataX1","Number of SPD Clusters in x, layer 1",
110                               9, 0.5, 9.5);
111     hdataX1->SetMarkerColor(2);
112     hdataX1->SetMarkerStyle(20);
113     hdataX1->SetMarkerSize(0.7);
114     TH1F *hbaSaX1  = new TH1F("hbaSaX1", "Ba/Sa", 11, 0.5, 11.5);
115     TH1F *hdubnaX1 = new TH1F("hdubnaX1", "Old Dubna", 7, 0.5, 7.5);
116     TH1F *hdataX2  = new TH1F("hdataX2","Number of SPD Clusters in x, layer 2",
117                               9, 0.5, 9.5);
118     hdataX2->SetMarkerColor(2);
119     hdataX2->SetMarkerStyle(20);
120     hdataX2->SetMarkerSize(0.7);
121     TH1F *hbaSaX2  = new TH1F("hbaSaX2", "Ba/Sa", 11, 0.5, 11.5);
122     TH1F *hdubnaX2 = new TH1F("hdubnaX2", "Old Dubna", 8, 0.5, 8.5);
123     TH1F *hdataZ1  = new TH1F("hdataZ1","Number of SPD Clusters in z, layer 1",
124                               4, 0.5, 4.5);
125     hdataZ1->SetMarkerColor(2);
126     hdataZ1->SetMarkerStyle(20);
127     hdataZ1->SetMarkerSize(0.7);
128     TH1F *hbaSaZ1  = new TH1F("hbaSaZ1", "Ba/Sa", 3, 0.5, 3.5);
129     TH1F *hdubnaZ1 = new TH1F("hdubnaZ1", "Old Dubna", 3, 0.5, 3.5);
130     TH1F *hdataZ2  = new TH1F("hdataZ2","Number of SPD Clusters in z, layer 2",
131                               4, 0.5, 4.5);
132     hdataZ2->SetMarkerColor(2);
133     hdataZ2->SetMarkerStyle(20);
134     hdataZ2->SetMarkerSize(0.7);
135     TH1F *hbaSaZ2  = new TH1F("hbaSaZ2", "Ba/Sa", 4, 0.5, 4.5);
136     TH1F *hdubnaZ2 = new TH1F("hdubnaZ2", "Old Dubna", 3, 0.5, 3.5);
137
138     Int_t j = 0; //j will be a counter for the loops to fill the histograms
139                  //with the values for clusters for the Data, the Ba/Sa model
140                  //and the old Dubna model
141     for(j=0; j<9; j++){
142         hdataX1->SetBinContent((j+1), (Double_t)dataX1[j]);
143         hdataX1->SetBinError((j+1), 0.00001);
144     }
145     for(j=0; j<11; j++) hbaSaX1->Fill((Float_t)j +0.5, baSaX1[j]);
146     for(j=0; j<7; j++)  hdubnaX1->Fill((Float_t)j +0.5, dubnaX1[j]);
147     for(j=0; j<9; j++){
148         hdataX2->SetBinContent((j+1), (Double_t)dataX2[j]);
149         hdataX2->SetBinError((j+1), 0.00001);
150     }
151     for(j=0; j<11; j++) hbaSaX2->Fill((Float_t)j +0.5, baSaX2[j]);
152     for(j=0; j<8; j++)  hdubnaX2->Fill((Float_t)j +0.5, dubnaX2[j]);
153     for(j=0; j<4; j++){
154         hdataZ1->SetBinContent((j+1), (Double_t)dataZ1[j]);
155         hdataZ1->SetBinError((j+1), 0.00001);
156     }
157     for(j=0; j<3; j++)  hbaSaZ1->Fill((Float_t)j +0.5, baSaZ1[j]);
158     for(j=0; j<3; j++)  hdubnaZ1->Fill((Float_t)j +0.5, dubnaZ1[j]);
159     for(j=0; j<4; j++){
160         hdataZ2->SetBinContent((j+1), (Double_t)dataZ2[j]);
161         hdataZ2->SetBinError((j+1), 0.00001);
162     }
163     for(j=0; j<4; j++)  hbaSaZ2->Fill((Float_t)j +0.5, baSaZ2[j]);
164     for(j=0; j<3; j++)  hdubnaZ2->Fill((Float_t)j +0.5, dubnaZ2[j]);
165
166     TTree *tc = 0;
167     TBranch *spdBranch = 0;
168     TClonesArray *spdClustcA = new TClonesArray("AliITSRawClusterSPD",1000);
169     AliITSRawClusterSPD *spdClust = 0;
170     char tn[20];
171     Int_t evnt,i,n,nc;
172     Float_t nclx = 0,nclz = 0;
173     for(evnt=0;evnt<nevents;evnt++){
174         sprintf(tn,"TreeC%d",evnt);
175         tc = (TTree*) file->Get(tn);
176         spdBranch = tc->GetBranch("ITSClustersSPD");
177         spdBranch->SetAddress(&spdClustcA);
178         n = (Int_t) tc->GetEntries();
179         for(i=0;i<n;i++){
180             tc->GetEvent(i);
181             nc = spdClustcA->GetEntries();
182             if(nc>240) nc = 240;
183             for(j=0;j<nc;j++){
184                 spdClust = (AliITSRawClusterSPD*) spdClustcA->At(j);
185                 nclx = spdClust->NclX();
186                 nclz = spdClust->NclZ();
187                 if(spdClust->Module()<80){
188                     hclx1->Fill(nclx-0.5);
189                     hclz1->Fill(nclz-0.5);
190                 }else if(spdClust->Module()<240){
191                     hclx2->Fill(nclx-0.5);
192                     hclz2->Fill(nclz-0.5);
193                 } //endif
194             } // end for j
195         } // end for i
196         spdClustcA->Clear();
197         delete spdBranch; spdBranch = 0;
198         delete tc; tc = 0;
199     } // end for evnt
200
201     //Begin Process of normalizing histograms
202     Double_t integral = (Double_t) hclx1->Integral();
203     if(integral>0.0) hclx1->Scale(1./integral);
204     integral = hclz1->Integral();
205     if(integral>0.0) hclz1->Scale(1./integral);
206     integral = hclx2->Integral();
207     if(integral>0.0) hclx2->Scale(1./integral);
208     integral = hclz2->Integral();
209     if(integral>0.0) hclz2->Scale(1./integral);
210
211     hdataX1->SetMinimum(0.000257);
212     hdataX1->SetMaximum(1.01);
213     hdataZ1->SetMinimum(0.000274);
214     hdataZ1->SetMaximum(1.01);
215     hdataX2->SetMinimum(0.000257);
216     hdataX2->SetMaximum(1.01);
217     hdataZ2->SetMinimum(0.000257);
218     hdataZ2->SetMaximum(1.01);
219
220     char    aDiffCoeff[50];            //Character array for sprintf
221     sprintf(aDiffCoeff,"SPD Clusters with Diffusion Coefficent=%f",diffCoeff);
222     
223     TCanvas *cSPDclusters = new TCanvas("cSPDclusters",aDiffCoeff,
224                                        400,10,600,700);
225     cSPDclusters->Divide(2, 2);
226
227     cSPDclusters->cd(1);
228     cSPDclusters->SetLogy();
229     hdataX1->Draw("e1p");
230     hclx1->Draw("same");
231     hbaSaX1->SetLineColor(4);
232     hbaSaX1->SetLineStyle(2);    
233     hbaSaX1->Draw("same");
234     hdubnaX1->SetLineColor(3);
235     hdubnaX1->SetLineStyle(4);
236     hdubnaX1->Draw("same");
237
238     cSPDclusters->cd(2);
239     cSPDclusters->SetLogy();
240     hdataZ1->Draw("e1p");
241     hclz1->Draw("same");
242     hbaSaZ2->SetLineColor(4);
243     hbaSaZ1->SetLineStyle(2);
244     hbaSaZ1->Draw("same");
245     hdubnaZ1->SetLineColor(3);
246     hdubnaZ1->SetLineStyle(4);
247     hdubnaZ1->Draw("same");
248
249     cSPDclusters->cd(3);
250     cSPDclusters->SetLogy();
251     hdataX2->Draw("e1p");
252     hclx2->Draw("same");
253     hbaSaX2->SetLineColor(4);
254     hbaSaX2->SetLineStyle(2);
255     hbaSaX2->Draw("same");
256     hdubnaX2->SetLineColor(3);
257     hdubnaX2->SetLineStyle(4);
258     hdubnaX2->Draw("same");
259
260     cSPDclusters->cd(4);
261     cSPDclusters->SetLogy();
262     hdataZ2->Draw("e1p");
263     hclz2->Draw("same");
264     hbaSaZ2->SetLineColor(4);
265     hbaSaZ2->SetLineStyle(2);
266     hbaSaZ2->Draw("same");
267     hdubnaZ2->SetLineColor(3);
268     hdubnaZ2->SetLineStyle(4);
269     hdubnaZ2->Draw("same");
270     cSPDclusters->cd();
271
272     if(gROOT->IsBatch()){
273         cSPDclusters->Print("SPDClusters.eps");
274     } // end if gROOT->IsBatch()
275
276 }