New macros for the offline shifts
[u/mrichter/AliRoot.git] / MUON / TestRecPoints.C
CommitLineData
cc727848 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id: TestRecPoints.C 23207 2007-12-20 09:59:20Z ivana $ */
17
18/// \ingroup macros
19/// \file TestRecPoints.C
20/// \brief This macro is to be used to check the trigger and tracking chamber behaviour
21/// through the RecPoints analysis
22///
23/// \author D. Stocco (Torino)
24
25#if !defined(__CINT__) || defined(__MAKECINT__)
26// ROOT
27#include "TString.h"
28#include "TH1.h"
29#include "TList.h"
30#include "TFile.h"
31#include "TDirectory.h"
32
33#include "TTree.h"
34#include "TClonesArray.h"
35
36// STEER
37#include "AliRun.h"
38#include "AliRunLoader.h"
39#include "AliLoader.h"
40#include "AliCDBManager.h"
41
42// tracker
43#include "AliGeomManager.h"
44
45// MUON
46#include "AliMpCDB.h"
47#include "AliMpDDL.h"
48#include "AliMpDDLStore.h"
49#include "AliMpConstants.h"
50
51// trigger
52#include "AliMUONVTriggerStore.h"
53#include "AliMUONDigitStoreV1.h"
54#include "AliMUONDigitMaker.h"
55#include "AliMUONLocalTrigger.h"
56
57// tracker
58#include "AliMUONClusterStoreV2.h"
59#include "AliMUONClusterFinderMLEM.h"
60#include "AliMUONPreClusterFinder.h"
61#include "AliMUONSimpleClusterServer.h"
62#include "AliMUONGeometryTransformer.h"
63#include "AliMUONVDigitStore.h"
64#include "AliMUONVCluster.h"
65
66#endif
67
68const Int_t kNtrigChambers = AliMpConstants::NofTriggerChambers();
69const Int_t kNcathodes = AliMpConstants::NofCathodes();
70
71const Int_t kNtrackChambers = AliMpConstants::NofTrackingChambers();
72
73enum {kOnlyTracker, kOnlyTrigger, kTrackTrig};
74
75Int_t GetPlane(Int_t ch, Int_t cath){return kNcathodes * ch + cath;}
76void ClusterSize(TList&, AliMUONVDigit*, Int_t&, Int_t);
77
78// Main Method
79void TestRecPoints(TString baseDir=".", TString outDir=".", Int_t whatToTest=kTrackTrig, Int_t runNumber=0, TString cdbStorage="local://$ALICE_ROOT")
80{
81 const Int_t kNplanes = kNtrigChambers * kNcathodes;
82 const Int_t kNslats = 18;
83
84 TList histoListTrig;
85
86 TH1F* trigStripMult[kNplanes];
87 TH1F* trigClusterSize[kNplanes];
88 TH1F* trigClusterMult[kNplanes];
89 TString histoName, histoTitle;
90
91 if(whatToTest!=kOnlyTracker){
92 for(Int_t ch=0; ch<kNtrigChambers; ch++){
93 for(Int_t cath=0; cath<kNcathodes; cath++){
94 Int_t iplane = GetPlane(ch, cath);
95 histoName = "stripMult";
96 histoName.Append(Form("Ch%iCath%i",ch,cath));
97 histoTitle = Form("Strip multiplicity Ch %i Cathode %i", ch, cath);
98 trigStripMult[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), kNslats, 0.-0.5, (Float_t)kNslats - 0.5);
99 trigStripMult[iplane]->SetXTitle("slat");
100 histoListTrig.Add(trigStripMult[iplane]);
101
102 histoName = "clusterSize";
103 histoName.Append(Form("Ch%iCath%i",ch,cath));
104 histoTitle = Form("Cluster size Ch %i Cathode %i", ch, cath);
105 trigClusterSize[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), 10, 0.-0.5, (Float_t)10 - 0.5);
106 trigClusterSize[iplane]->SetXTitle("cluster size");
107 histoListTrig.Add(trigClusterSize[iplane]);
108
109 histoName = "clusterMult";
110 histoName.Append(Form("Ch%iCath%i",ch,cath));
111 histoTitle = Form("Cluster multiplicity Ch %i Cathode %i", ch, cath);
112 trigClusterMult[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), kNslats, 0.-0.5, (Float_t)kNslats - 0.5);
113 trigClusterMult[iplane]->SetXTitle("slat");
114 histoListTrig.Add(trigClusterMult[iplane]);
115 }
116 }
117 }
118
119 const Int_t kMaxDetElem = 26;
120 TList histoListTrack;
121 TH1F* trackClusterSize[kNtrackChambers];
122 TH1F* trackClusterMult[kNtrackChambers];
123 TH1F* trackClusterYield[kNtrackChambers];
124
125 if(whatToTest!=kOnlyTrigger){
126 for(Int_t ch=0; ch<kNtrackChambers; ch++){
127 histoName = "clusterSize";
128 histoName.Append(Form("Ch%i",ch));
129 histoTitle = Form("Cluster size Ch %i", ch);
130 trackClusterSize[ch] = new TH1F(histoName.Data(), histoTitle.Data(), kMaxDetElem, 0.-0.5, (Float_t)kMaxDetElem - 0.5);
131 trackClusterSize[ch]->SetXTitle("cluster size");
132 histoListTrack.Add(trackClusterSize[ch]);
133
134 histoName = "clusterMult";
135 histoName.Append(Form("Ch%i",ch));
136 histoTitle = Form("Cluster multiplicity vs. Charge Ch %i", ch);
137 trackClusterMult[ch] = new TH1F(histoName.Data(), histoTitle.Data(), 200, 0, 15000.);
138 trackClusterMult[ch]->SetXTitle("Charge (pC)");
139 histoListTrack.Add(trackClusterMult[ch]);
140
141 histoName = "clusterYield";
142 histoName.Append(Form("Ch%i",ch));
143 histoTitle = Form("Cluster yield vs. DetElem Ch %i", ch);
144 trackClusterYield[ch] = new TH1F(histoName.Data(), histoTitle.Data(), kMaxDetElem, 0.-0.5, (Float_t)kMaxDetElem - 0.5);
145 trackClusterYield[ch]->SetXTitle("Detector element #");
146 histoListTrack.Add(trackClusterYield[ch]);
147 }
148 }
149
150
151 // Creating Run Loader and opening RecPoints
152 TString filename = baseDir.Data();
153 filename.Append("/galice.root");
154 AliRunLoader * RunLoader = AliRunLoader::Open(filename.Data(),"MUONLoader","UPDATE");
155 if (RunLoader ==0x0) {
156 printf(">>> Error : Error Opening %s file \n",filename.Data());
157 return;
158 }
159
160 // Loading MUON subsystem
161 AliLoader* MUONLoader = RunLoader->GetDetectorLoader("MUON");
162 if(whatToTest!=kOnlyTracker) MUONLoader->LoadRecPoints("READ");
163 if(whatToTest!=kOnlyTrigger) MUONLoader->LoadDigits("READ");
164
165 TTree *treeR = 0x0, *treeD = 0x0;
166 AliMUONVTriggerStore* trigStore = 0x0;
167 AliMUONLocalTrigger* locTrg = 0x0;
168
169 AliMUONVDigitStore* digitStoreTrack = 0x0;
170 AliMUONVCluster* cluster = 0x0;
171
172 // Load segmentation
173 AliCDBManager::Instance()->SetDefaultStorage(cdbStorage.Data());
174 AliCDBManager::Instance()->SetRun(runNumber);
175
176 AliMpCDB::LoadDDLStore();
177
178 AliMUONGeometryTransformer* transformer = 0x0;
179
180
181 AliMUONClusterStoreV2 clusterStore;
182 AliMUONVClusterFinder* clusterFinder = 0x0;
183
184 AliMUONSimpleClusterServer* clusterServer = 0x0;
185
186 Int_t firstChamber(0);
187 Int_t lastChamber(9);
188
189 if(whatToTest!=kOnlyTrigger){
190 // Import TGeo geometry
191 TString geoFilename = baseDir.Data();
192 geoFilename.Append("/geometry.root");
193 if ( ! AliGeomManager::GetGeometry() ) {
194 AliGeomManager::LoadGeometry(geoFilename);
195 if (! AliGeomManager::GetGeometry() ) {
196 printf(">>> Error : getting geometry from file %s failed\n", geoFilename.Data());
197 return;
198 }
199 }
200 transformer = new AliMUONGeometryTransformer();
201 // Load geometry data
202 transformer->LoadGeometryData();
203 clusterFinder = new AliMUONClusterFinderMLEM(kFALSE,new AliMUONPreClusterFinder);
204 clusterServer = new AliMUONSimpleClusterServer(*clusterFinder,*transformer);
205 }
206
207 AliMUONDigitStoreV1 digitStore;
208 AliMUONVDigit* mDigit;
209
210 Int_t clusterSize;
211
212 AliMUONDigitMaker digitMaker;
213 TList digitsList[kNplanes];
214
215 Int_t nevents = RunLoader->GetNumberOfEvents();
216 Int_t analysisFrac = nevents/10+1;
217
218 printf("\nNumber of events = %i\n\n",nevents);
219
220 for(Int_t ievent=0; ievent<nevents; ievent++){
221 if(ievent%analysisFrac==0) printf("Analysing event = %i\n",ievent);
222 RunLoader->GetEvent(ievent);
223 if(whatToTest!=kOnlyTracker){
224 digitStore.Clear();
225 treeR = MUONLoader->TreeR();
226 trigStore = AliMUONVTriggerStore::Create(*treeR);
227
228 if ( trigStore == 0x0 ) continue;
229 trigStore->Clear();
230 trigStore->Connect(*treeR);
231 treeR->GetEvent(0);
232
233 TIter nextLocal(trigStore->CreateLocalIterator());
234 while ( (locTrg = static_cast<AliMUONLocalTrigger*>( nextLocal() )) != NULL )
235 {
236 TArrayS xyPattern[2];
237 locTrg->GetXPattern(xyPattern[0]);
238 locTrg->GetYPattern(xyPattern[1]);
239
240 Int_t nBoard = locTrg->LoCircuit();
241 digitMaker.TriggerDigits(nBoard, xyPattern, digitStore);
242 }
243
244 TIter next(digitStore.CreateIterator());
245
246 while ( ( mDigit = static_cast<AliMUONVDigit*>(next()) ) )
247 {
248 Int_t detElemId = mDigit->DetElemId();
249 Int_t ch = detElemId/100 - 11;
250 Int_t cathode = mDigit->Cathode();
251 Int_t slat = detElemId%100;
252 Int_t iplane = GetPlane(ch, cathode);
253 trigStripMult[iplane]->Fill(slat);
254 digitsList[iplane].Add(mDigit);
255 } // loop on digits
256
257 for(Int_t iplane=0; iplane<kNplanes; iplane++){
258 while(digitsList[iplane].GetEntries()){
259 clusterSize=1;
260 mDigit = (AliMUONVDigit*)digitsList[iplane].At(0);
261 digitsList[iplane].Remove(mDigit);
262 ClusterSize(digitsList[iplane],mDigit,clusterSize,1);
263 //if(clusterSize>1) printf("Cluster size = %i\n\n",clusterSize);
264 trigClusterSize[iplane]->Fill(clusterSize);
265
266 Int_t detElemId = mDigit->DetElemId();
267 Int_t slat = detElemId%100;
268 trigClusterMult[iplane]->Fill(slat);
269 } // loop o sorted digits
270 } // loop on planes
271 } // trigger part
272
273 if(whatToTest!=kOnlyTrigger){
274 clusterStore.Clear();
275 treeD = MUONLoader->TreeD();
276 digitStoreTrack = AliMUONVDigitStore::Create(*treeD);
277
278 if ( digitStoreTrack == 0x0 ) continue;
279 digitStoreTrack->Clear();
280 digitStoreTrack->Connect(*treeD);
281 treeD->GetEvent(0);
282 clusterServer->UseDigitStore(*digitStoreTrack);
283
284 for (Int_t ch = firstChamber; ch <= lastChamber; ++ch ){
285 clusterServer->Clusterize(ch, clusterStore, AliMpArea());
286 }
287
288 TIter next(clusterStore.CreateIterator());
289 while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) )
290 {
291 Float_t charge = cluster->GetCharge();
292 if(charge==0.) continue;
293 Int_t ch = cluster->GetChamberId();
294 Int_t npads = cluster->GetNDigits();
295 Int_t detElemId = cluster->GetDetElemId()%100;
296 //printf("charge = %f pads = %i detElemId = %i\n",charge,npads,detElemId);
297
298 trackClusterSize[ch]->Fill(npads);
299 trackClusterMult[ch]->Fill(charge);
300 trackClusterYield[ch]->Fill(detElemId);
301 } // loop on clusters
302 } // tracker part
303 } // loop on events
304
305 MUONLoader->UnloadRecPoints();
306 MUONLoader->UnloadDigits();
307 RunLoader->UnloadAll();
308 delete RunLoader;
309 TString outFileName = outDir.Data();
310 outFileName.Append("/outTestRecPoints.root");
311 TFile* outFile = new TFile(outFileName.Data(), "RECREATE");
312 TDirectory* dir = 0x0;
313 if(whatToTest!=kOnlyTracker){
314 outFile->cd();
315 dir = outFile->mkdir("Trigger");
316 dir->cd();
317 histoListTrig.Write();
318 }
319 if(whatToTest!=kOnlyTrigger){
320 outFile->cd();
321 dir = outFile->mkdir("Tracker");
322 dir->cd();
323 histoListTrack.Write();
324 }
325 outFile->Close();
326 printf("\nSee results in %s\n",outFileName.Data());
327}
328
329void ClusterSize(TList& list, AliMUONVDigit* refDigit, Int_t& clusterSize, Int_t depthLevel)
330{
331 AliMUONVDigit* mDigit = 0x0;
332 for(Int_t idigit=0; idigit<list.GetEntries(); idigit++){
333 mDigit = (AliMUONVDigit*) list.At(idigit);
334 if(mDigit->DetElemId() != refDigit->DetElemId()) continue;
335 Int_t diffX = TMath::Abs(mDigit->PadX() - refDigit->PadX());
336 Int_t diffY = TMath::Abs(mDigit->PadY() - refDigit->PadY());
337 if(diffX>1) continue;
338 if(diffY>1) continue;
339 if(diffX*diffY != 0) continue;
340 clusterSize++;
341 list.Remove(mDigit);
342 //printf("DetElemId = %i Level = %i Ref. (%2i,%2i) Found (%2i,%2i)\n",mDigit->DetElemId(),depthLevel,refDigit->PadX(),refDigit->PadY(),mDigit->PadX(),mDigit->PadY());
343 ClusterSize(list, mDigit, clusterSize, depthLevel+1);
344 Int_t val = idigit + depthLevel - clusterSize;
345 //printf("Level = %i Current digit = %i\t",depthLevel,idigit);
346 idigit = TMath::Max(-1,val);
347 //printf("New digit = %i\n",idigit);
348 }
349}