Changes related to the extraction of the V0 finder into a separate class (A. Dainese...
[u/mrichter/AliRoot.git] / ITS / AliITSAnalizeSPDHits.C
1 void AliITSAnalizeSPDHits(TString hfn="galice.root",Int_t mod=-1,
2                      Int_t evnt=-1){
3   // Macro to analize hits in the SPD to chatch posible problems
4
5   // Dynamically link some shared libs
6   if (gClassTable->GetID("AliRun") < 0) {
7     gROOT->LoadMacro("loadlibs.C");
8     loadlibs();
9   } else {
10     if(gAlice){
11       delete AliRunLoader::Instance();
12       delete gAlice;
13       gAlice=0;
14     }
15   }
16   gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSstandard.C");
17   // Set OCDB if needed
18   AliCDBManager* man = AliCDBManager::Instance();
19   if (!man->IsDefaultStorageSet()) {
20     printf("Setting a local default storage and run number 0\n");
21     man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
22     man->SetRun(0);
23   }else {
24     printf("Using deafult storage \n");
25   }
26   // retrives geometry 
27   TString geof(gSystem->DirName(hfn));
28   geof += "/geometry.root";
29   TGeoManager::Import(geof.Data());
30   if (!gGeoManager) {
31     cout<<"geometry not found\n";
32     return -1;
33   }
34
35
36   AliRunLoader *rl = AccessFile(hfn); // Set up to read in Data
37   Int_t retval = rl->LoadHeader();
38   if (retval){
39     cerr<<"AliITSPrintHits.C : LoadHeader returned error"<<endl;
40     return;
41   }
42
43   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
44
45   if(!ITSloader){
46     cerr<<"AliITSPrintHits.C :  ITS loader not found"<<endl;
47     return;
48   }
49
50   ITSloader->LoadHits("read");
51   AliITS *ITS = (AliITS*)gAlice->GetDetector("ITS");
52   if(!ITS || ITS==0){
53     cout << "Error: no ITS found. Aborting"<<endl;
54     return;
55   } // end if !ITS
56   //cout << ITS << endl;
57   AliITSInitGeometry *initgeom = new AliITSInitGeometry(ITS->GetMajorVersion(),ITS->GetMinorVersion());
58   //cout << initgeom << endl;
59   AliITSgeom *geom = initgeom->CreateAliITSgeom();
60   if(!geom){
61     cout << "Error: not AliITSgeom object found."<<endl;
62     return;
63   } // end if geom
64   ITSloader->SetITSgeom(geom);
65
66   Int_t evNumber1 = 0;
67   Int_t evNumber2 = AliRunLoader::GetNumberOfEvents();
68   if(evnt>=0){
69     evNumber1 = evnt;
70     evNumber2 = evnt+1;
71   } // end if evnt>=0
72   Int_t mod1 = 0;
73   Int_t mod2 = geom->GetIndexMax();
74   if(mod>=0){
75     mod1 = mod;
76     mod2 = mod+1;
77   } // end if mod>=0
78   AliITShit *hp = 0;
79
80   Int_t i;
81   Double_t zbin[161];
82   for(i=0;i<161;i++) zbin[i] = 425.0E-4;
83   zbin[32] = zbin[64] = zbin[96] = zbin[128] = 625.0E-4;
84   zbin[0] = -3.536;
85   for(i=1;i<161;i++) zbin[i] += zbin[i-1];
86   TH1I *phiDistribution1 = new TH1I("phiDistribution1",
87                    "Phi Distribution for the inner most layer of the SPD",
88                                     5027,0.0,360.0);
89   TH2I *detCoverage1 = new TH2I("detCoverage1",
90               "The SPD layer 1 cell coverage summed over the first ladder per sector",
91                                 160,(Double_t *)zbin,256,-0.64,+0.64);
92   TH2I *detCoverage1a = new TH2I("detCoverage1a",
93               "The SPD layer 1 cell coverage sum over the second ladder per sector",
94                                  160,(Double_t *)zbin,256,-0.64,+0.64);
95   TH2I *detThicknessZ1 = new TH2I("detThicknessZ1",
96                     "Hit local y distribution as a function of z",
97                                   160,(Double_t*)zbin,200,-100.E-4,+100.E-4);
98   TH2I *detThicknessX1 = new TH2I("detThicknessX1",
99                     "Hit local y distribution as a function of x",
100                                   256,-0.64,+0.64,200,-100.E-4,+100.E-4);
101   TH1I *phiDistribution2 = new TH1I("phiDistribution2",
102                    "Phi Distribution for the outer most layer of the SPD",
103                                     8800,0.0,360.0);
104   TH2I *detCoverage2 = new TH2I("detCoverage2",
105               "The SPD layer 2 cell coverage summed over all but the forth ladder per sector",
106                                 160,(Double_t *)zbin,256,-0.64,+0.64);
107   TH2I *detCoverage2a = new TH2I("detCoverage2a",
108               "The SPD layer 2 cell coverage summed over all forth ladder per sector",
109                                  160,(Double_t *)zbin,256,-0.64,+0.64);
110   TH2I *detThicknessZ2 = new TH2I("detThicknessZ2",
111                     "Hit local y distribution as a function of z",
112                                   160,(Double_t*)zbin,200,-100.E-4,+100.E-4);
113   TH2I *detThicknessX2 = new TH2I("detThicknessX2",
114                     "Hit local y distribution as a function of x",
115                                   256,-0.64,+0.64,200,-100.E-4,+100.E-4);
116   TH2I *xySpace = new TH2I("xySpace",
117                            "ALICE global coorinate location of hits in X,y",
118                            100,-8.0,+8.0,100,-8.0,+8.0);
119   TH2I *zrSpace = new TH2I("zrSpace",
120                            "ALICE global coorinate location of hits in z,r",
121                            100,-15.0,+15.0,100,0.0,+8.0);
122   Double_t xg,yg,zg,xl,yl,zl,phi;
123   Int_t nmodules,size=-1;
124   Int_t event,m,i,i2,hit,trk,lay,lad,det;
125   for(event = evNumber1; event < evNumber2; event++){
126     //cout<<"Processing event "<<event<<endl;
127     rl->GetEvent(event);
128     ITS->InitModules(size,nmodules);
129     ITS->FillModules(event,0,-1," "," ");
130     for(m=mod1;m<mod2;m++)if(geom->GetModuleType(m)==(AliITSDetector)(0)){
131       i2 = (ITS->GetModule(m))->GetNhits();
132       //cout <<  "Event=" << event << " module=" << m <<
133       //  " Number of Hits=" << i2 <<endl;
134       for(i=0;i<i2;i++){
135         trk = (ITS->GetModule(m))->GetHitTrackIndex(i);
136         hit = (ITS->GetModule(m))->GetHitHitIndex(i);
137         hp  = (ITS->GetModule(m))->GetHit(i);
138         hp->GetPositionG(xg,yg,zg);
139         hp->GetPositionL(xl,yl,zl);
140         geom->GetModuleId(m,lay,lad,det);
141         phi = TMath::ATan2(yg,xg)*TMath::RadToDeg();
142         if(phi<0.0) phi+=360.0;
143         switch(lay){
144         case 1:
145           phiDistribution1->Fill(phi,1.0);
146           if(lad%2==1) detCoverage1->Fill(zl,xl,1.0);
147           if(lad%2==0) detCoverage1a->Fill(zl,xl,1.0);
148           detThicknessZ1->Fill(zl,yl,1.0);
149           detThicknessX1->Fill(xl,yl,1.0);
150           break;
151         case 2:
152           phiDistribution2->Fill(phi,1.0);
153           if(lad%4!=0) detCoverage2->Fill(zl,xl,1.0);
154           if(lad%4==0) detCoverage2a->Fill(zl,xl,1.0);
155           detThicknessZ2->Fill(zl,yl,1.0);
156           detThicknessX2->Fill(xl,yl,1.0);
157         break;
158         default:
159         } // end switch
160       } // end for i
161     } // end for m
162     ITS->ClearModules();
163   } // end for event
164   //
165   cout << "Creating pad"<<endl;
166   TVirtualPad *pad=0;
167   TCanvas *c0 = new TCanvas("c0","SPD Hits Tests",262,50,700,534);
168   c0->Range(0,0,1,1);
169   c0->Divide(3,2);
170   //
171   pad = c0->cd(1);
172   phiDistribution2->Draw();
173   pad = c0->cd(4);
174   phiDistribution1->Draw();
175   TPad *pad2 = c0->cd(2);
176   delete pad2;
177   c0->cd(0);
178   TPad *pad21 = new TPad("pad21","CoverageZ",0.34,0.875,0.65,1.0);
179   pad21->Draw();
180   pad21->cd();
181   pad21->Range(0,0,1,1);
182   pad21->Modified();
183   detCoverage2->Draw();
184   c0->Update();
185   c0->cd(0);
186   TPad *pad22 = new TPad("pad22","CoverageZa",0.34,0.75,0.65,0.875);
187   pad22->Draw();
188   pad22->cd();
189   pad22->Range(0,0,1,1);
190   pad22->Modified();
191   detCoverage2a->Draw();
192   c0->Update();
193   c0->cd(0);
194   TPad *pad23 = new TPad("pad23","ThicknessZ",0.34,0.625,0.65,0.75);
195   pad23->Draw();
196   pad23->cd();
197   pad23->Range(0,0,1,1);
198   pad23->Modified();
199   detThicknessZ2->Draw();
200   c0->Update();
201   c0->cd(0);
202   TPad *pad24 = new TPad("pad24","ThicknessX",0.34,0.50,0.65,0.625);
203   pad24->Draw();
204   pad24->cd();
205   pad24->Range(0,0,1,1);
206   pad24->Modified();
207   detThicknessX2->Draw();
208   c0->Update();
209   //
210   c0->cd(0);
211   TPad *pad5 = c0->cd(5);
212   delete pad5;
213   c0->cd(0);
214   TPad *pad51 = new TPad("pad51","CoverageZ",0.34,0.375,0.65,0.50);
215   pad51->Draw();
216   pad51->cd();
217   pad51->Range(0,0,1,1);
218   pad51->Modified();
219   detCoverage1->Draw();
220   c0->Update();
221   c0->cd(0);
222   TPad *pad52 = new TPad("pad52","CoverageZa",0.34,0.25,0.65,0.375);
223   pad52->Draw();
224   pad52->cd();
225   pad52->Range(0,0,1,1);
226   pad52->Modified();
227   detCoverage1a->Draw();
228   c0->Update();
229   c0->cd(0);
230   TPad *pad53 = new TPad("pad53","ThicknessZ",0.34,0.125,0.65,0.25);
231   pad53->Draw();
232   pad53->cd();
233   pad53->Range(0,0,1,1);
234   pad53->Modified();
235   detThicknessZ1->Draw();
236   c0->Update();
237   c0->cd(0);
238   TPad *pad54 = new TPad("pad54","ThicknessX",0.34,0.00,0.65,0.125);
239   pad54->Draw();
240   pad54->cd();
241   pad54->Range(0,0,1,1);
242   pad54->Modified();
243   detThicknessX1->Draw();
244   c0->Update();
245   //
246   TList *tvTreeList = new TList;
247   TTree *tvTree = (TTree *) gROOT->FindObject("TreeH");
248   tvTreeList->Add(tvTree);
249   pad = c0->cd(3);
250   tvTree->Draw("ITS.fY:ITS.fX","ITS.fModule<240","", 100000, 0);
251   pad = c0->cd(6);
252   tvTree->Draw("ITS.fX*ITS.fX+ITS.fY*ITS.fY:ITS.fZ","ITS.fModule<240","", 
253                  100000, 0);
254
255 }