Make the Scan method public
[u/mrichter/AliRoot.git] / MUON / TestRecPoints.C
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 "AliMpConstants.h"
47
48 // trigger
49 #include "AliMUONVTriggerStore.h"
50 #include "AliMUONDigitStoreV1.h"
51 #include "AliMUONDigitMaker.h"
52 #include "AliMUONLocalTrigger.h"
53
54 // tracker
55 #include "AliMUONClusterStoreV2.h"
56 #include "AliMUONClusterFinderMLEM.h"
57 #include "AliMUONPreClusterFinder.h"
58 #include "AliMUONSimpleClusterServer.h"
59 #include "AliMUONGeometryTransformer.h"
60 #include "AliMUONVDigitStore.h"
61 #include "AliMUONDigitStoreV2R.h"
62 #include "AliMUONVCluster.h"
63 #include "AliMUONRecoParam.h"
64 #include "AliMUONCDB.h"
65
66 #endif
67
68 const Int_t kNtrigChambers = AliMpConstants::NofTriggerChambers();
69 const Int_t kNcathodes = AliMpConstants::NofCathodes();
70
71 const Int_t kNtrackChambers = AliMpConstants::NofTrackingChambers();
72
73 enum {kOnlyTracker, kOnlyTrigger, kTrackTrig};
74
75 Int_t GetPlane(Int_t ch, Int_t cath){return kNcathodes * ch + cath;}
76 void ClusterSize(TList&, AliMUONVDigit*, Int_t&, Int_t);
77
78 // Main Method
79 void TestRecPoints(TString baseDir=".", TString outDir=".", Float_t adcCut = 10., Int_t whatToTest=kTrackTrig, Int_t runNumber=0, TString cdbStorage="local://$ALICE_ROOT/OCDB")
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 (ADC)");
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   AliMUONDigitStoreV2R digitStoreTrackCut;
171   AliMUONVCluster* cluster = 0x0;
172   
173   // Load mapping
174   AliCDBManager::Instance()->SetDefaultStorage(cdbStorage.Data());
175   AliCDBManager::Instance()->SetRun(runNumber);
176   if (!AliMUONCDB::LoadMapping()) return;
177   
178   AliMUONGeometryTransformer* transformer = 0x0;
179   
180   AliMUONRecoParam* recoParam = 0x0;
181   
182   AliMUONClusterStoreV2 clusterStore;
183   AliMUONVClusterFinder* clusterFinder = 0x0;
184   
185   AliMUONSimpleClusterServer* clusterServer = 0x0;
186   
187   Int_t firstChamber(0);
188   Int_t lastChamber(9);
189   
190   if(whatToTest!=kOnlyTrigger){
191     // Import TGeo geometry
192     TString geoFilename = baseDir.Data();
193     geoFilename.Append("/geometry.root");
194     if ( ! AliGeomManager::GetGeometry() ) {
195       AliGeomManager::LoadGeometry(geoFilename);
196       if (! AliGeomManager::GetGeometry() ) {
197         printf(">>> Error : getting geometry from file %s failed\n", geoFilename.Data());
198         return;
199       }
200     }
201     transformer = new AliMUONGeometryTransformer();
202     // Load geometry data
203     transformer->LoadGeometryData();
204     // Load reconstruction parameters
205     recoParam = AliMUONCDB::LoadRecoParam();
206     if (!recoParam) return;
207     recoParam->Print("FULL");
208     clusterFinder = new AliMUONClusterFinderMLEM(kFALSE,new AliMUONPreClusterFinder);
209     clusterServer = new AliMUONSimpleClusterServer(clusterFinder,*transformer);
210   }
211   
212   AliMUONDigitStoreV1 digitStore;
213   AliMUONVDigit* mDigit;
214   
215   Int_t clusterSize;
216   
217   AliMUONDigitMaker digitMaker;
218   TList digitsList[kNplanes];
219   
220   Int_t nevents = RunLoader->GetNumberOfEvents();
221   Int_t analysisFrac = nevents/10+1;
222   
223   printf("\nNumber of events = %i\n\n",nevents);
224   
225   for(Int_t ievent=0; ievent<nevents; ievent++){
226     if(ievent%analysisFrac==0) printf("Analysing event = %i\n",ievent);
227     RunLoader->GetEvent(ievent);
228     if(whatToTest!=kOnlyTracker){
229       digitStore.Clear();
230       treeR = MUONLoader->TreeR();
231       trigStore = AliMUONVTriggerStore::Create(*treeR);
232       
233       if ( trigStore == 0x0 ) continue;
234       trigStore->Clear();
235       trigStore->Connect(*treeR);
236       treeR->GetEvent(0);
237       
238       TIter nextLocal(trigStore->CreateLocalIterator());
239       while ( (locTrg = static_cast<AliMUONLocalTrigger*>( nextLocal() )) != NULL )
240       {
241         TArrayS xyPattern[2];
242         locTrg->GetXPattern(xyPattern[0]);
243         locTrg->GetYPattern(xyPattern[1]);
244         
245         Int_t nBoard = locTrg->LoCircuit();
246         digitMaker.TriggerDigits(nBoard, xyPattern, digitStore);
247       }
248       
249       TIter next(digitStore.CreateIterator());
250       
251       while ( ( mDigit = static_cast<AliMUONVDigit*>(next()) ) )
252       {
253         Int_t detElemId = mDigit->DetElemId();
254         Int_t ch = detElemId/100 - 11;
255         Int_t cathode = mDigit->Cathode();
256         Int_t slat = detElemId%100;
257         Int_t iplane = GetPlane(ch, cathode);
258         trigStripMult[iplane]->Fill(slat);
259         digitsList[iplane].Add(mDigit);
260       } // loop on digits
261       
262       for(Int_t iplane=0; iplane<kNplanes; iplane++){
263         while(digitsList[iplane].GetEntries()){
264           clusterSize=1;
265           mDigit = (AliMUONVDigit*)digitsList[iplane].At(0);
266           digitsList[iplane].Remove(mDigit);
267           ClusterSize(digitsList[iplane],mDigit,clusterSize,1);
268           //if(clusterSize>1) printf("Cluster size = %i\n\n",clusterSize);
269           trigClusterSize[iplane]->Fill(clusterSize);
270           
271           Int_t detElemId = mDigit->DetElemId();
272           Int_t slat = detElemId%100;
273           trigClusterMult[iplane]->Fill(slat);
274         } // loop o sorted digits
275       } // loop on planes
276     } // trigger part
277     
278     if(whatToTest!=kOnlyTrigger){
279       clusterStore.Clear();
280       treeD = MUONLoader->TreeD();
281       digitStoreTrack = AliMUONVDigitStore::Create(*treeD);
282       
283       if ( digitStoreTrack == 0x0 ) continue;
284       digitStoreTrack->Clear();
285       digitStoreTrackCut.Clear();
286       digitStoreTrack->Connect(*treeD);
287       treeD->GetEvent(0);
288       
289       // Cut low charge channels (pedestal subtraction)
290       TIter nextDigitTrack(digitStoreTrack->CreateIterator());
291       while ( ( mDigit = static_cast<AliMUONVDigit*>(nextDigitTrack()) ) )
292       {
293         //printf("Digit class = %s",mDigit->ClassName());
294         Float_t charge = mDigit->Charge();
295         if(charge<adcCut) continue;
296         digitStoreTrackCut.Add(*mDigit, AliMUONVDigitStore::kDeny);
297       } // loop on digits
298       
299       TIter nextDigitTrackCut(digitStoreTrackCut.CreateIterator());
300       clusterServer->UseDigits(nextDigitTrackCut,&digitStoreTrackCut);
301       
302       for (Int_t ch = firstChamber; ch <= lastChamber; ++ch ){
303         clusterServer->Clusterize(ch, clusterStore, AliMpArea(),recoParam);
304       }
305       
306       TIter next(clusterStore.CreateIterator());
307       while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) )
308       {
309         Float_t charge = cluster->GetCharge();
310         if(charge==0.) continue;
311         Int_t ch = cluster->GetChamberId();
312         Int_t npads = cluster->GetNDigits();
313         Int_t detElemId = cluster->GetDetElemId()%100;
314         //printf("charge = %f   pads = %i   detElemId = %i\n",charge,npads,detElemId);
315         
316         trackClusterSize[ch]->Fill(npads);
317         trackClusterMult[ch]->Fill(charge);
318         trackClusterYield[ch]->Fill(detElemId);      
319       } // loop on clusters
320     } // tracker part
321   } // loop on events
322   
323   MUONLoader->UnloadRecPoints();
324   MUONLoader->UnloadDigits();
325   RunLoader->UnloadAll();
326   delete RunLoader;
327   TString outFileName = outDir.Data();
328   outFileName.Append("/outTestRecPoints.root");
329   TFile* outFile = new TFile(outFileName.Data(), "RECREATE");
330   TDirectory* dir = 0x0;
331   if(whatToTest!=kOnlyTracker){
332     outFile->cd();
333     dir = outFile->mkdir("Trigger");
334     dir->cd();
335     histoListTrig.Write();
336   }
337   if(whatToTest!=kOnlyTrigger){
338     outFile->cd();
339     dir = outFile->mkdir("Tracker");
340     dir->cd();
341     histoListTrack.Write();
342   }
343   outFile->Close();
344   printf("\nSee results in %s\n",outFileName.Data());
345 }
346
347 void ClusterSize(TList& list, AliMUONVDigit* refDigit, Int_t& clusterSize, Int_t depthLevel)
348 {
349   AliMUONVDigit* mDigit = 0x0;
350   for(Int_t idigit=0; idigit<list.GetEntries(); idigit++){
351     mDigit = (AliMUONVDigit*) list.At(idigit);
352     if(mDigit->DetElemId() != refDigit->DetElemId()) continue;
353     Int_t diffX = TMath::Abs(mDigit->PadX() - refDigit->PadX());
354     Int_t diffY = TMath::Abs(mDigit->PadY() - refDigit->PadY());
355     if(diffX>1) continue;
356     if(diffY>1) continue;
357     if(diffX*diffY != 0) continue;
358     clusterSize++;
359     list.Remove(mDigit);
360     //printf("DetElemId = %i   Level = %i   Ref. (%2i,%2i)    Found (%2i,%2i)\n",mDigit->DetElemId(),depthLevel,refDigit->PadX(),refDigit->PadY(),mDigit->PadX(),mDigit->PadY());
361     ClusterSize(list, mDigit, clusterSize, depthLevel+1);
362     Int_t val = idigit + depthLevel - clusterSize;
363     //printf("Level = %i   Current digit = %i\t",depthLevel,idigit);
364     idigit = TMath::Max(-1,val);
365     //printf("New digit = %i\n",idigit);
366   }
367 }