Adding new macros in loadmacros.C
[u/mrichter/AliRoot.git] / MUON / TestRecPoints.C
CommitLineData
cc727848 1/**************************************************************************
6f6a3f3b 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**************************************************************************/
cc727848 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"
35be7ed7 41#include "AliCDBEntry.h"
42#include "AliCDBPath.h"
cc727848 43
44// tracker
45#include "AliGeomManager.h"
46
47// MUON
48#include "AliMpCDB.h"
49#include "AliMpDDL.h"
50#include "AliMpDDLStore.h"
51#include "AliMpConstants.h"
52
53// trigger
54#include "AliMUONVTriggerStore.h"
55#include "AliMUONDigitStoreV1.h"
56#include "AliMUONDigitMaker.h"
57#include "AliMUONLocalTrigger.h"
58
59// tracker
60#include "AliMUONClusterStoreV2.h"
61#include "AliMUONClusterFinderMLEM.h"
62#include "AliMUONPreClusterFinder.h"
63#include "AliMUONSimpleClusterServer.h"
64#include "AliMUONGeometryTransformer.h"
65#include "AliMUONVDigitStore.h"
6f6a3f3b 66#include "AliMUONDigitStoreV2R.h"
cc727848 67#include "AliMUONVCluster.h"
35be7ed7 68#include "AliMUONRecoParam.h"
cc727848 69
70#endif
71
72const Int_t kNtrigChambers = AliMpConstants::NofTriggerChambers();
73const Int_t kNcathodes = AliMpConstants::NofCathodes();
74
75const Int_t kNtrackChambers = AliMpConstants::NofTrackingChambers();
76
77enum {kOnlyTracker, kOnlyTrigger, kTrackTrig};
78
79Int_t GetPlane(Int_t ch, Int_t cath){return kNcathodes * ch + cath;}
80void ClusterSize(TList&, AliMUONVDigit*, Int_t&, Int_t);
81
82// Main Method
162637e4 83void TestRecPoints(TString baseDir=".", TString outDir=".", Float_t adcCut = 10., Int_t whatToTest=kTrackTrig, Int_t runNumber=0, TString cdbStorage="local://$ALICE_ROOT/OCDB")
cc727848 84{
85 const Int_t kNplanes = kNtrigChambers * kNcathodes;
86 const Int_t kNslats = 18;
6f6a3f3b 87
cc727848 88 TList histoListTrig;
6f6a3f3b 89
cc727848 90 TH1F* trigStripMult[kNplanes];
91 TH1F* trigClusterSize[kNplanes];
92 TH1F* trigClusterMult[kNplanes];
93 TString histoName, histoTitle;
6f6a3f3b 94
cc727848 95 if(whatToTest!=kOnlyTracker){
96 for(Int_t ch=0; ch<kNtrigChambers; ch++){
97 for(Int_t cath=0; cath<kNcathodes; cath++){
6f6a3f3b 98 Int_t iplane = GetPlane(ch, cath);
99 histoName = "stripMult";
100 histoName.Append(Form("Ch%iCath%i",ch,cath));
101 histoTitle = Form("Strip multiplicity Ch %i Cathode %i", ch, cath);
102 trigStripMult[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), kNslats, 0.-0.5, (Float_t)kNslats - 0.5);
103 trigStripMult[iplane]->SetXTitle("slat");
104 histoListTrig.Add(trigStripMult[iplane]);
105
106 histoName = "clusterSize";
107 histoName.Append(Form("Ch%iCath%i",ch,cath));
108 histoTitle = Form("Cluster size Ch %i Cathode %i", ch, cath);
109 trigClusterSize[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), 10, 0.-0.5, (Float_t)10 - 0.5);
110 trigClusterSize[iplane]->SetXTitle("cluster size");
111 histoListTrig.Add(trigClusterSize[iplane]);
112
113 histoName = "clusterMult";
114 histoName.Append(Form("Ch%iCath%i",ch,cath));
115 histoTitle = Form("Cluster multiplicity Ch %i Cathode %i", ch, cath);
116 trigClusterMult[iplane] = new TH1F(histoName.Data(), histoTitle.Data(), kNslats, 0.-0.5, (Float_t)kNslats - 0.5);
117 trigClusterMult[iplane]->SetXTitle("slat");
118 histoListTrig.Add(trigClusterMult[iplane]);
cc727848 119 }
120 }
121 }
6f6a3f3b 122
cc727848 123 const Int_t kMaxDetElem = 26;
124 TList histoListTrack;
125 TH1F* trackClusterSize[kNtrackChambers];
126 TH1F* trackClusterMult[kNtrackChambers];
127 TH1F* trackClusterYield[kNtrackChambers];
6f6a3f3b 128
cc727848 129 if(whatToTest!=kOnlyTrigger){
130 for(Int_t ch=0; ch<kNtrackChambers; ch++){
131 histoName = "clusterSize";
132 histoName.Append(Form("Ch%i",ch));
133 histoTitle = Form("Cluster size Ch %i", ch);
134 trackClusterSize[ch] = new TH1F(histoName.Data(), histoTitle.Data(), kMaxDetElem, 0.-0.5, (Float_t)kMaxDetElem - 0.5);
135 trackClusterSize[ch]->SetXTitle("cluster size");
136 histoListTrack.Add(trackClusterSize[ch]);
6f6a3f3b 137
cc727848 138 histoName = "clusterMult";
139 histoName.Append(Form("Ch%i",ch));
140 histoTitle = Form("Cluster multiplicity vs. Charge Ch %i", ch);
141 trackClusterMult[ch] = new TH1F(histoName.Data(), histoTitle.Data(), 200, 0, 15000.);
6f6a3f3b 142 trackClusterMult[ch]->SetXTitle("Charge (ADC)");
cc727848 143 histoListTrack.Add(trackClusterMult[ch]);
6f6a3f3b 144
cc727848 145 histoName = "clusterYield";
146 histoName.Append(Form("Ch%i",ch));
147 histoTitle = Form("Cluster yield vs. DetElem Ch %i", ch);
148 trackClusterYield[ch] = new TH1F(histoName.Data(), histoTitle.Data(), kMaxDetElem, 0.-0.5, (Float_t)kMaxDetElem - 0.5);
149 trackClusterYield[ch]->SetXTitle("Detector element #");
150 histoListTrack.Add(trackClusterYield[ch]);
151 }
152 }
153
6f6a3f3b 154
cc727848 155 // Creating Run Loader and opening RecPoints
156 TString filename = baseDir.Data();
157 filename.Append("/galice.root");
158 AliRunLoader * RunLoader = AliRunLoader::Open(filename.Data(),"MUONLoader","UPDATE");
159 if (RunLoader ==0x0) {
160 printf(">>> Error : Error Opening %s file \n",filename.Data());
161 return;
162 }
6f6a3f3b 163
cc727848 164 // Loading MUON subsystem
165 AliLoader* MUONLoader = RunLoader->GetDetectorLoader("MUON");
166 if(whatToTest!=kOnlyTracker) MUONLoader->LoadRecPoints("READ");
167 if(whatToTest!=kOnlyTrigger) MUONLoader->LoadDigits("READ");
6f6a3f3b 168
cc727848 169 TTree *treeR = 0x0, *treeD = 0x0;
170 AliMUONVTriggerStore* trigStore = 0x0;
171 AliMUONLocalTrigger* locTrg = 0x0;
6f6a3f3b 172
cc727848 173 AliMUONVDigitStore* digitStoreTrack = 0x0;
6f6a3f3b 174 AliMUONDigitStoreV2R digitStoreTrackCut;
cc727848 175 AliMUONVCluster* cluster = 0x0;
6f6a3f3b 176
cc727848 177 // Load segmentation
178 AliCDBManager::Instance()->SetDefaultStorage(cdbStorage.Data());
179 AliCDBManager::Instance()->SetRun(runNumber);
6f6a3f3b 180
cc727848 181 AliMpCDB::LoadDDLStore();
6f6a3f3b 182
cc727848 183 AliMUONGeometryTransformer* transformer = 0x0;
6f6a3f3b 184
35be7ed7 185 AliMUONRecoParam* recoParam = 0x0;
6f6a3f3b 186
cc727848 187 AliMUONClusterStoreV2 clusterStore;
188 AliMUONVClusterFinder* clusterFinder = 0x0;
189
190 AliMUONSimpleClusterServer* clusterServer = 0x0;
6f6a3f3b 191
cc727848 192 Int_t firstChamber(0);
193 Int_t lastChamber(9);
6f6a3f3b 194
cc727848 195 if(whatToTest!=kOnlyTrigger){
196 // Import TGeo geometry
197 TString geoFilename = baseDir.Data();
198 geoFilename.Append("/geometry.root");
199 if ( ! AliGeomManager::GetGeometry() ) {
200 AliGeomManager::LoadGeometry(geoFilename);
201 if (! AliGeomManager::GetGeometry() ) {
6f6a3f3b 202 printf(">>> Error : getting geometry from file %s failed\n", geoFilename.Data());
203 return;
cc727848 204 }
205 }
206 transformer = new AliMUONGeometryTransformer();
207 // Load geometry data
208 transformer->LoadGeometryData();
35be7ed7 209 // Load reconstruction parameters
210 AliCDBPath path("MUON","Calib","RecoParam");
211 AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
212 if(entry) {
213 recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
214 entry->SetOwner(0);
215 AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
216 }
217 if (!recoParam) {
218 printf("Couldn't find RecoParam object in OCDB: create default one");
219 recoParam = AliMUONRecoParam::GetLowFluxParam();
220 }
221 recoParam->Print("FULL");
cc727848 222 clusterFinder = new AliMUONClusterFinderMLEM(kFALSE,new AliMUONPreClusterFinder);
0f211038 223 clusterServer = new AliMUONSimpleClusterServer(clusterFinder,*transformer);
cc727848 224 }
6f6a3f3b 225
cc727848 226 AliMUONDigitStoreV1 digitStore;
227 AliMUONVDigit* mDigit;
6f6a3f3b 228
cc727848 229 Int_t clusterSize;
6f6a3f3b 230
cc727848 231 AliMUONDigitMaker digitMaker;
232 TList digitsList[kNplanes];
6f6a3f3b 233
cc727848 234 Int_t nevents = RunLoader->GetNumberOfEvents();
235 Int_t analysisFrac = nevents/10+1;
6f6a3f3b 236
cc727848 237 printf("\nNumber of events = %i\n\n",nevents);
6f6a3f3b 238
cc727848 239 for(Int_t ievent=0; ievent<nevents; ievent++){
240 if(ievent%analysisFrac==0) printf("Analysing event = %i\n",ievent);
241 RunLoader->GetEvent(ievent);
242 if(whatToTest!=kOnlyTracker){
243 digitStore.Clear();
244 treeR = MUONLoader->TreeR();
245 trigStore = AliMUONVTriggerStore::Create(*treeR);
6f6a3f3b 246
cc727848 247 if ( trigStore == 0x0 ) continue;
248 trigStore->Clear();
249 trigStore->Connect(*treeR);
250 treeR->GetEvent(0);
6f6a3f3b 251
cc727848 252 TIter nextLocal(trigStore->CreateLocalIterator());
253 while ( (locTrg = static_cast<AliMUONLocalTrigger*>( nextLocal() )) != NULL )
254 {
6f6a3f3b 255 TArrayS xyPattern[2];
256 locTrg->GetXPattern(xyPattern[0]);
257 locTrg->GetYPattern(xyPattern[1]);
258
259 Int_t nBoard = locTrg->LoCircuit();
260 digitMaker.TriggerDigits(nBoard, xyPattern, digitStore);
cc727848 261 }
6f6a3f3b 262
cc727848 263 TIter next(digitStore.CreateIterator());
6f6a3f3b 264
cc727848 265 while ( ( mDigit = static_cast<AliMUONVDigit*>(next()) ) )
266 {
6f6a3f3b 267 Int_t detElemId = mDigit->DetElemId();
268 Int_t ch = detElemId/100 - 11;
269 Int_t cathode = mDigit->Cathode();
270 Int_t slat = detElemId%100;
271 Int_t iplane = GetPlane(ch, cathode);
272 trigStripMult[iplane]->Fill(slat);
273 digitsList[iplane].Add(mDigit);
cc727848 274 } // loop on digits
6f6a3f3b 275
cc727848 276 for(Int_t iplane=0; iplane<kNplanes; iplane++){
6f6a3f3b 277 while(digitsList[iplane].GetEntries()){
278 clusterSize=1;
279 mDigit = (AliMUONVDigit*)digitsList[iplane].At(0);
280 digitsList[iplane].Remove(mDigit);
281 ClusterSize(digitsList[iplane],mDigit,clusterSize,1);
282 //if(clusterSize>1) printf("Cluster size = %i\n\n",clusterSize);
283 trigClusterSize[iplane]->Fill(clusterSize);
284
285 Int_t detElemId = mDigit->DetElemId();
286 Int_t slat = detElemId%100;
287 trigClusterMult[iplane]->Fill(slat);
288 } // loop o sorted digits
cc727848 289 } // loop on planes
290 } // trigger part
6f6a3f3b 291
cc727848 292 if(whatToTest!=kOnlyTrigger){
293 clusterStore.Clear();
294 treeD = MUONLoader->TreeD();
295 digitStoreTrack = AliMUONVDigitStore::Create(*treeD);
6f6a3f3b 296
cc727848 297 if ( digitStoreTrack == 0x0 ) continue;
298 digitStoreTrack->Clear();
6f6a3f3b 299 digitStoreTrackCut.Clear();
cc727848 300 digitStoreTrack->Connect(*treeD);
301 treeD->GetEvent(0);
6f6a3f3b 302
303 // Cut low charge channels (pedestal subtraction)
304 TIter nextDigitTrack(digitStoreTrack->CreateIterator());
305 while ( ( mDigit = static_cast<AliMUONVDigit*>(nextDigitTrack()) ) )
306 {
307 //printf("Digit class = %s",mDigit->ClassName());
308 Float_t charge = mDigit->Charge();
309 if(charge<adcCut) continue;
310 digitStoreTrackCut.Add(*mDigit, AliMUONVDigitStore::kDeny);
311 } // loop on digits
312
0f211038 313 TIter nextDigitTrackCut(digitStoreTrackCut.CreateIterator());
2e2d0c44 314 clusterServer->UseDigits(nextDigitTrackCut,&digitStoreTrackCut);
6f6a3f3b 315
cc727848 316 for (Int_t ch = firstChamber; ch <= lastChamber; ++ch ){
35be7ed7 317 clusterServer->Clusterize(ch, clusterStore, AliMpArea(),recoParam);
cc727848 318 }
6f6a3f3b 319
cc727848 320 TIter next(clusterStore.CreateIterator());
321 while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) )
322 {
6f6a3f3b 323 Float_t charge = cluster->GetCharge();
324 if(charge==0.) continue;
325 Int_t ch = cluster->GetChamberId();
326 Int_t npads = cluster->GetNDigits();
327 Int_t detElemId = cluster->GetDetElemId()%100;
328 //printf("charge = %f pads = %i detElemId = %i\n",charge,npads,detElemId);
329
330 trackClusterSize[ch]->Fill(npads);
331 trackClusterMult[ch]->Fill(charge);
332 trackClusterYield[ch]->Fill(detElemId);
cc727848 333 } // loop on clusters
334 } // tracker part
335 } // loop on events
6f6a3f3b 336
cc727848 337 MUONLoader->UnloadRecPoints();
338 MUONLoader->UnloadDigits();
339 RunLoader->UnloadAll();
340 delete RunLoader;
341 TString outFileName = outDir.Data();
342 outFileName.Append("/outTestRecPoints.root");
343 TFile* outFile = new TFile(outFileName.Data(), "RECREATE");
344 TDirectory* dir = 0x0;
345 if(whatToTest!=kOnlyTracker){
346 outFile->cd();
347 dir = outFile->mkdir("Trigger");
348 dir->cd();
349 histoListTrig.Write();
350 }
351 if(whatToTest!=kOnlyTrigger){
352 outFile->cd();
353 dir = outFile->mkdir("Tracker");
354 dir->cd();
355 histoListTrack.Write();
356 }
357 outFile->Close();
358 printf("\nSee results in %s\n",outFileName.Data());
359}
360
361void ClusterSize(TList& list, AliMUONVDigit* refDigit, Int_t& clusterSize, Int_t depthLevel)
362{
363 AliMUONVDigit* mDigit = 0x0;
364 for(Int_t idigit=0; idigit<list.GetEntries(); idigit++){
365 mDigit = (AliMUONVDigit*) list.At(idigit);
366 if(mDigit->DetElemId() != refDigit->DetElemId()) continue;
367 Int_t diffX = TMath::Abs(mDigit->PadX() - refDigit->PadX());
368 Int_t diffY = TMath::Abs(mDigit->PadY() - refDigit->PadY());
369 if(diffX>1) continue;
370 if(diffY>1) continue;
371 if(diffX*diffY != 0) continue;
372 clusterSize++;
373 list.Remove(mDigit);
374 //printf("DetElemId = %i Level = %i Ref. (%2i,%2i) Found (%2i,%2i)\n",mDigit->DetElemId(),depthLevel,refDigit->PadX(),refDigit->PadY(),mDigit->PadX(),mDigit->PadY());
375 ClusterSize(list, mDigit, clusterSize, depthLevel+1);
376 Int_t val = idigit + depthLevel - clusterSize;
377 //printf("Level = %i Current digit = %i\t",depthLevel,idigit);
378 idigit = TMath::Max(-1,val);
379 //printf("New digit = %i\n",idigit);
380 }
381}