Updates in order to enable the '2D' PID for the TRD developed by Daniel Lohner.
[u/mrichter/AliRoot.git] / ANALYSIS / TenderSupplies / AliTRDTenderSupply.cxx
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
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 // TRD tender: reapply pid on the fly                                        //
20 //                                                                           //
21 ///////////////////////////////////////////////////////////////////////////////
22
23 #include <TChain.h>
24 #include <TDirectory.h>
25 #include <TFile.h>
26 #include <TList.h>
27 #include <TObjString.h>
28 #include <TTree.h>
29 #include <TString.h>
30 #include <TVectorD.h>
31
32 #include <AliCDBEntry.h>
33 #include <AliCDBId.h>
34 #include <AliCDBManager.h>
35 #include <AliOADBContainer.h>
36 #include <AliTRDCalDet.h>
37
38 #include <AliLog.h>
39 #include <TTree.h>
40 #include <TChain.h>
41 #include <AliGeomManager.h>
42 #include <AliPID.h>
43 #include <AliVEvent.h>
44 #include <AliESDEvent.h>
45 #include <AliESDpid.h>
46 #include <AliESDtrack.h>
47 #include <AliESDInputHandler.h>
48 #include <AliAnalysisManager.h>
49 #include <AliTrackerBase.h>
50 #include <AliTRDPIDResponseObject.h>
51 #include <AliTRDPIDResponse.h>
52 #include "AliTRDCalChamberStatus.h"
53 #include <AliTender.h>
54
55 #include "AliTRDTenderSupply.h"
56
57 ClassImp(AliTRDTenderSupply)
58
59 //_____________________________________________________
60 AliTRDTenderSupply::AliTRDTenderSupply() :
61   AliTenderSupply(),
62   fESD(NULL),
63   fESDpid(NULL),
64   fChamberGainOld(NULL),
65   fChamberGainNew(NULL),
66   fChamberVdriftOld(NULL),
67   fChamberVdriftNew(NULL),
68   fRunByRunCorrection(NULL),
69   fPIDmethod(kNNpid),
70   fNormalizationFactor(1.),
71   fPthreshold(0.8),
72   fNBadChambers(0),
73   fGeoFile(NULL),
74   fGainCorrection(kTRUE),
75   fLoadReferences(kFALSE),
76   fLoadReferencesFromCDB(kFALSE),
77   fLoadDeadChambers(kFALSE),
78   fHasReferences(kFALSE),
79   fHasNewCalibration(kTRUE),
80   fDebugMode(kFALSE),
81   fNameRunByRunCorrection()
82 {
83   //
84   // default ctor
85   //
86   memset(fSlicesForPID, 0, sizeof(UInt_t) * 2);
87   memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers);
88 }
89
90 //_____________________________________________________
91 AliTRDTenderSupply::AliTRDTenderSupply(const char *name, const AliTender *tender) :
92   AliTenderSupply(name,tender),
93   fESD(NULL),
94   fESDpid(NULL),
95   fChamberGainOld(NULL),
96   fChamberGainNew(NULL),
97   fChamberVdriftOld(NULL),
98   fChamberVdriftNew(NULL),
99   fRunByRunCorrection(NULL),
100   fPIDmethod(kNNpid),
101   fNormalizationFactor(1.),
102   fPthreshold(0.8),
103   fNBadChambers(0),
104   fGeoFile(NULL),
105   fGainCorrection(kTRUE),
106   fLoadReferences(kFALSE),
107   fLoadReferencesFromCDB(kFALSE),
108   fLoadDeadChambers(kFALSE),
109   fHasReferences(kFALSE),
110   fHasNewCalibration(kTRUE),
111   fDebugMode(kFALSE),
112   fNameRunByRunCorrection()
113 {
114   //
115   // named ctor
116   //
117   memset(fSlicesForPID, 0, sizeof(UInt_t) * 2);
118   memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers);
119 }
120
121 //_____________________________________________________
122 AliTRDTenderSupply::~AliTRDTenderSupply()
123 {
124   //
125   // dtor
126   //
127 }
128
129 //_____________________________________________________
130 void AliTRDTenderSupply::Init()
131 {
132   //
133   // Initialise TRD tender
134   //
135
136   AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
137   
138   // 1DLQ PID implemented in the AliESD object
139   fESDpid=fTender->GetESDhandler()->GetESDpid();
140   if (!fESDpid) {
141     fESDpid=new AliESDpid;
142     fTender->GetESDhandler()->SetESDpid(fESDpid);
143   }
144   // Load References
145   if(fLoadReferences && !fLoadReferencesFromCDB) LoadReferences();
146   //fESDpid->GetTRDResponse().SetGainNormalisationFactor(fNormalizationFactor);
147   fESDpid->SetTRDslicesForPID(fSlicesForPID[0], fSlicesForPID[1]);
148
149   if(fNameRunByRunCorrection.Length()) LoadRunByRunCorrection(fNameRunByRunCorrection.Data());
150
151   // Set Normalisation Factors
152   if(mgr->GetMCtruthEventHandler()){
153     // Assume MC
154     //fESDpid->GetTRDResponse().SetGainNormalisationFactor(1.);
155     SwitchOffGainCorrection();
156   }
157   else{
158     // Assume Data
159     //if(fPIDmethod == kNNpid) fPidRecal->SetGainScaleFactor(1.14);
160     //fESDpid->GetTRDResponse().SetGainNormalisationFactor(1.14);
161     SwitchOnGainCorrection();
162   }
163 }
164
165 //_____________________________________________________
166 void AliTRDTenderSupply::ProcessEvent()
167 {
168   //
169   // Reapply pid information
170   //
171   if (fTender->RunChanged()){
172     AliDebug(0, Form("AliTPCTenderSupply::ProcessEvent - Run Changed (%d)\n",fTender->GetRun()));
173     if (fGainCorrection) SetChamberGain();
174     if(fLoadReferences && !fHasReferences) LoadReferences();
175     if(fLoadDeadChambers) LoadDeadChambersFromCDB();
176     // Load Geometry
177     if(AliGeomManager::GetGeometry()){
178       AliInfo("Geometry already loaded by other tenders");
179     } else {
180       if(fGeoFile) AliInfo(Form("Load geometry from file %s\n", fGeoFile));
181       else AliInfo("Load Geometry from OCDB\n");
182       AliGeomManager::LoadGeometry(fGeoFile);
183     }
184   }
185
186
187   fESD = fTender->GetEvent();
188   if (!fESD) return;
189   Int_t ntracks=fESD->GetNumberOfTracks();
190    
191   //
192   // recalculate PID probabilities
193   //
194   Int_t detectors[kNPlanes];
195   for(Int_t itrack = 0; itrack < ntracks; itrack++){
196     for(Int_t idet = 0; idet < 5; idet++) detectors[idet] = -1;
197     AliESDtrack *track=fESD->GetTrack(itrack);
198     // Recalculate likelihoods
199     if(!(track->GetStatus() & AliESDtrack::kTRDout)) continue;
200     AliDebug(2, Form("TRD track found, gain correction: %s, Number of bad chambers: %d\n", fGainCorrection ? "Yes" : "No", fNBadChambers));
201     if(GetTRDchamberID(track, detectors)){
202       if(fGainCorrection && fHasNewCalibration) ApplyGainCorrection(track, detectors);
203       if(fNBadChambers) MaskChambers(track, detectors);
204     }
205     if(fRunByRunCorrection) ApplyRunByRunCorrection(track);
206     if(fNormalizationFactor != 1.){
207       //printf("Gain Factor: %f\n", fNormalizationFactor);
208       // Renormalize charge
209       Double_t qslice = -1;
210       for(Int_t ily = 0; ily < 6; ily++){
211         for(Int_t is = 0; is < track->GetNumberOfTRDslices(); is++){
212           qslice = track->GetTRDslice(ily, is);
213           //printf("Doing layer %d slice %d, value %f\n", ily, is, qslice);
214           if(qslice >0){
215             qslice *= fNormalizationFactor;
216             //printf("qslice new: %f\n", qslice);
217             track->SetTRDslice(qslice, ily, is);
218           }
219         }
220       }
221     }
222     switch(fPIDmethod){
223       case kNNpid:
224         break;
225       case k1DLQpid:
226         fESDpid->MakeTRDPID(fESD->GetTrack(itrack));
227         break;
228       default:
229         AliError("PID Method not implemented (yet)");
230     }
231   }
232 }
233
234 //_____________________________________________________
235 void AliTRDTenderSupply::LoadDeadChambersFromCDB(){
236   //
237   // Load Dead Chambers from the OCDB
238   //
239   AliDebug(1, "Loading Dead Chambers from the OCDB");
240   AliCDBEntry *en = fTender->GetCDBManager()->Get("TRD/Calib/ChamberStatus",fTender->GetRun());
241   if(!en){
242    AliError("Dead Chambers not in OCDB");
243    return;
244   }
245   en->GetId().Print();
246
247   AliTRDCalChamberStatus* chamberStatus = 0;
248   if(en){
249     chamberStatus = (AliTRDCalChamberStatus*)en->GetObject();
250     if(!chamberStatus) AliError("List with the dead chambers not found");
251     for(Int_t ichamber = 0; ichamber < 540; ichamber++) {
252       if(!chamberStatus->IsGood(ichamber)){
253         //printf("Chamber not installed %d\n",ichamber);
254         AddBadChamber(ichamber);
255       }
256     }
257   }
258 }
259
260 //_____________________________________________________
261 void AliTRDTenderSupply::LoadReferences(){
262   //
263   // Load Reference from the OCDB/OADB into the PID Response
264   //
265   if(fLoadReferencesFromCDB){
266     AliDebug(1, "Loading Reference Distributions from the OCDB");
267     AliCDBEntry *en = fTender->GetCDBManager()->Get("TRD/Calib/PIDLQ1D");
268     if(!en){
269       AliError("References for 1D Likelihood Method not in OCDB");
270       return;
271     }
272     en->GetId().Print();
273     TObjArray *arr = dynamic_cast<TObjArray *>(en->GetObject());
274     if(!arr) AliError("List with the references not found");
275   
276     // Get new references
277     TIter refs(arr);
278     TObject *o = NULL;
279     AliTRDPIDResponseObject *ref = NULL;
280     while((o = refs())){
281       if(!TString(o->IsA()->GetName()).CompareTo("AliTRDPIDResponseObject")){
282         ref = dynamic_cast<AliTRDPIDResponseObject *>(o);
283         break;
284       }
285     }
286     if(ref){
287       fESDpid->GetTRDResponse().SetPIDResponseObject(ref);
288       fHasReferences = kTRUE;
289       AliDebug(1, "Reference distributions loaded into the PID Response");
290     } else {
291       AliError("References not found");
292     }
293   } else {
294     // Backward compatibility mode
295     AliInfo("Loading Reference Distributions from ROOT file");
296     fESDpid->GetTRDResponse().Load("$TRAIN_ROOT/util/tender/LQ1dRef_v3.root");
297     fHasReferences = kTRUE;
298   }
299 }
300
301 //_____________________________________________________
302 void AliTRDTenderSupply::SetChamberGain(){
303   //
304   // Load Chamber Gain factors into the Tender supply
305   //
306   
307   //find previous entry from the UserInfo
308   TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree();
309   if (!tree) {
310     fHasNewCalibration = kFALSE;
311     AliError("Tree not found in ESDhandler");
312     return;
313   }
314          
315   TList *userInfo=(TList*)tree->GetUserInfo();
316   if (!userInfo) {
317     fHasNewCalibration = kFALSE;
318     AliError("No UserInfo found in tree");
319     return;
320   }
321
322   TList *cdbList=(TList*)userInfo->FindObject("cdbList");
323   if (!cdbList) {
324     fHasNewCalibration = kFALSE;
325     AliError("No cdbList found in UserInfo");
326     if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
327     return;
328   }
329   fHasNewCalibration = kTRUE;
330         
331   TIter nextCDB(cdbList);
332   TObjString *os=0x0;
333   while ( (os=(TObjString*)nextCDB()) ){
334     if(os->GetString().Contains("TRD/Calib/ChamberGainFactor")){
335       // Get Old gain calibration
336       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
337            
338       AliCDBEntry *entry=fTender->GetCDBManager()->Get(id->GetPath(), id->GetFirstRun(), id->GetVersion());
339       if (!entry) {
340         AliError("No previous gain calibration entry found");
341         return;
342       }
343
344       fChamberGainOld = dynamic_cast<AliTRDCalDet *>(entry->GetObject());
345            
346       AliDebug(1, Form("Used old Gain entry: %s\n",entry->GetId().ToString().Data()));
347     } else if(os->GetString().Contains("TRD/Calib/ChamberVdrift")){
348       // Get Old drift velocity calibration
349       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
350            
351       AliCDBEntry *entry=fTender->GetCDBManager()->Get(id->GetPath(), id->GetFirstRun(), id->GetVersion());
352       if (!entry) {
353         AliError("No previous drift velocity calibration entry found");
354         return;
355       }
356
357       fChamberVdriftOld = dynamic_cast<AliTRDCalDet *>(entry->GetObject());
358            
359       AliDebug(1, Form("Used old Drift velocity entry: %s\n",entry->GetId().ToString().Data()));
360     
361     }
362   }
363
364   // Get Latest Gain Calib Object
365   AliCDBEntry *entryNew=fTender->GetCDBManager()->Get("TRD/Calib/ChamberGainFactor",fTender->GetRun());
366   if (entryNew) {
367     AliDebug(1, Form("Used new Gain entry: %s\n",entryNew->GetId().ToString().Data()));
368     fChamberGainNew = dynamic_cast<AliTRDCalDet *>(entryNew->GetObject());
369   } else
370     AliError("No new gain calibration entry found");
371   
372   // Also get the latest Drift Velocity calibration object
373   entryNew=fTender->GetCDBManager()->Get("TRD/Calib/ChamberVdrift",fTender->GetRun());
374   if (entryNew) {
375     AliDebug(1, Form("Used new Drift velocity entry: %s\n",entryNew->GetId().ToString().Data()));
376     fChamberVdriftNew = dynamic_cast<AliTRDCalDet *>(entryNew->GetObject());
377   } else
378     AliError("No new drift velocity calibration entry found");
379
380   if(!fChamberGainNew || !fChamberVdriftNew){
381     AliError("No recent calibration found");
382     fHasNewCalibration = kFALSE;
383   }
384 }
385
386 //_____________________________________________________
387 void AliTRDTenderSupply::LoadRunByRunCorrection(const char *filename){
388   //
389   // Define run by run gain correction for the charge
390   // 
391
392   TDirectory *bkp = gDirectory;
393   TFile *in = TFile::Open(filename);
394   bkp->cd();
395   fRunByRunCorrection = dynamic_cast<AliOADBContainer *>(in->Get("TRDchargeCorrection"));
396   delete in;
397   if(fRunByRunCorrection )
398     AliDebug(2, Form("OADB Container has %d runs\n", fRunByRunCorrection->GetNumberOfEntries()));
399   /* Temporarily out due to a bug in AliOADBContainer
400   fRunByRunCorrection = new AliOADBContainer("TRDchargeCorrection");
401   Int_t status = fRunByRunCorrection->InitFromFile(filename, "TRDchargeCorrection");
402   if(!status) AliDebug(1, Form("Run-dependend gain correction factors loaded from OADB file %s", filename));
403   else{
404     AliDebug(1, "Failed Loading Run-dependend gain correction factors");
405     delete fRunByRunCorrection;
406     fRunByRunCorrection = NULL;
407   }
408   */
409 }
410
411 //_____________________________________________________
412 Bool_t AliTRDTenderSupply::IsBadChamber(Int_t chamberID){
413   //
414   // Check if the chamber id is in the list of bad chambers
415   //
416   Bool_t isBad = kFALSE;
417   for(UInt_t icam = 0; icam < fNBadChambers; icam++)
418     if(fBadChamberID[icam] == chamberID){
419             isBad = kTRUE;
420             //printf("cross checking: %i \n",chamberID);
421       break;
422     }
423   return isBad;
424 }
425
426 //_____________________________________________________
427 void AliTRDTenderSupply::ApplyGainCorrection(AliESDtrack * track, const Int_t * const chamberID){
428   //
429   // Apply new gain factors to the track
430   //
431   if(!fChamberGainNew || !fChamberGainOld){
432     AliError("Cannot apply gain correction.");
433     return;
434   }
435  
436   if(!(track->GetStatus() & AliESDtrack::kTRDout)) return;
437   Bool_t applyCorrectionVdrift = kFALSE;
438   if(fChamberVdriftOld && fChamberVdriftNew) applyCorrectionVdrift = kTRUE;
439
440   for(Int_t iplane = 0; iplane < kNPlanes; iplane++){
441     if(chamberID[iplane] < 0) continue;
442     if(IsBadChamber(chamberID[iplane])) continue; // Don't apply gain correction for chambers which are in the list of bad chambers
443
444     // Take old and new gain factor and make ratio
445     Double_t facOld = fChamberGainOld->GetValue(chamberID[iplane]);
446     Double_t facNew = fChamberGainNew->GetValue(chamberID[iplane]); 
447     Double_t correction = facNew/facOld;
448     if(applyCorrectionVdrift){
449       // apply also correction for drift velocity calibration
450       Double_t vDriftOld = fChamberVdriftOld->GetValue(chamberID[iplane]);
451       Double_t vDriftNew = fChamberVdriftNew->GetValue(chamberID[iplane]);
452       correction *= vDriftNew/vDriftOld;
453     }
454     AliDebug(2, Form("Applying correction factor %f\n", correction));
455     for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
456       Double_t qslice = track->GetTRDslice(iplane, islice);
457       if(qslice <= 0.) continue; 
458       track->SetTRDslice(qslice / correction, iplane, islice);
459     }
460   }
461 }
462
463 //_____________________________________________________
464 void AliTRDTenderSupply::ApplyRunByRunCorrection(AliESDtrack *const track) {
465   // 
466   // Equalize charge distribution by applying run-by-run correction (multiplicative)
467   //
468
469   TVectorD *corrfactor = dynamic_cast<TVectorD *>(fRunByRunCorrection->GetObject(fTender->GetRun()));
470   if(!corrfactor) {
471     AliDebug(2, "Couldn't derive gain correction factor from OADB");
472     return;
473   }
474   else AliDebug(2, Form("Gain factor from OADB %f", (*corrfactor)[0]));
475   Double_t slice = 0;
476   for(Int_t ily = 0; ily < kNPlanes; ily++){
477     for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
478       slice = track->GetTRDslice(ily, islice);
479       if(slice < 0.001) continue;   // do not modify slices which are 0 or negative
480       slice *= (*corrfactor)[0];
481       track->SetTRDslice(slice, ily, islice);
482     }
483   }
484 }
485
486 //_____________________________________________________
487 void AliTRDTenderSupply::MaskChambers(AliESDtrack *const track, const Int_t * const chamberID){
488   //
489   // Mask out chambers which are in the list of bad chambers
490   // Set chamber signal to 0 and reduce the number of tracklets used for PID
491   //
492   AliDebug(2, "Masking bad chambers for TRD track");
493   Int_t nTrackletsPID = 0, nslice = 0, nTracklets = track->GetTRDntracklets();
494   Bool_t badChamber = kFALSE;
495   //Int_t nbad = 0 ;      // Number of bad chambers which contain also a signal
496   //Int_t nsliceBad = 0;  // Number of slices in tracklet in a bad chamber
497   for(Int_t iplane = 0; iplane < kNPlanes; iplane++){
498     badChamber = kFALSE;
499     nslice = 0; //nsliceBad = 0;
500     if(IsBadChamber(chamberID[iplane])) badChamber = kTRUE;
501     for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
502       if(badChamber){
503         //if(track->GetTRDslice(iplane, islice)) nsliceBad++;
504         track->SetTRDslice(-1, iplane, islice);
505       } else if(track->GetTRDslice(iplane, islice) > 0.001) nslice++;
506     }
507     //if(nsliceBad) nbad++;
508     if(nslice > 0) nTrackletsPID++;
509   }
510   //if(nbad) track->SetTRDncls(track->GetTRDncls() - 20 * nbad);      // subtract mean number of clusters per tracklet for bad tracklets 
511   // Use nTrackletsPID to indicate the number of tracklets from good
512   // chambers so they are used for the PID
513   track->SetTRDntracklets(nTrackletsPID | (nTracklets << 3));
514 }
515
516 //_____________________________________________________
517 Bool_t AliTRDTenderSupply::GetTRDchamberID(AliESDtrack * const track, Int_t *detectors) {
518   //
519   // Calculate TRD chamber ID
520   //
521   Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
522   if(p < fPthreshold) return kFALSE; // Apply low momentum cutoff
523
524   Double_t xLayer[kNPlanes] = {300.2, 312.8, 325.4, 338., 350.6, 363.2};
525   Double_t etamin[kNStacks] = {0.536, 0.157, -0.145, -0.527,-0.851};
526   Double_t etamax[kNStacks] = {0.851, 0.527, 0.145, -0.157,-0.536};
527   for(Int_t ily = 0; ily < kNPlanes; ily++) detectors[ily] = -1;
528
529   const AliExternalTrackParam *trueparam = NULL;
530   if(track->GetOuterParam()) trueparam = track->GetOuterParam();
531   else if(track->GetTPCInnerParam()) trueparam = track->GetTPCInnerParam();
532   else if(track->GetInnerParam()) trueparam = track->GetInnerParam();
533   if(!trueparam){
534     AliDebug(2, "No Track Params");
535     return kFALSE;
536   }
537
538   AliExternalTrackParam workparam(*trueparam); // Do calculation on working Copy
539   Double_t pos[3];
540   Int_t nDet = 0;
541   for(Int_t ily = 0; ily < kNPlanes; ily++){
542     if(!AliTrackerBase::PropagateTrackToBxByBz(&workparam, xLayer[ily], 0.139, 100)){   // Assuming the pion mass
543       AliDebug(2, "Propagation failed");
544       break;
545     }
546     workparam.GetXYZ(pos);
547     Double_t trackAlpha = TMath::ATan2(pos[1], pos[0]);
548     if(fDebugMode){
549       // Compare to simple propagation without magnetic field
550       AliExternalTrackParam workparam1(*trueparam); // Do calculation on working Copy
551       Double_t pos1[3];
552       if(!workparam1.PropagateTo(xLayer[ily], fESD->GetMagneticField())) {
553         AliDebug(2, "Propagation failed");
554         break;
555       }
556       workparam1.GetXYZ(pos1);
557       Double_t trackAlpha1 = TMath::ATan2(pos1[1], pos1[0]);
558       AliDebug(2, Form("Alpha: Old %f, New %f, diff %f", trackAlpha1, trackAlpha, trackAlpha-trackAlpha1));
559     }
560     if(trackAlpha < 0) trackAlpha = 2 * TMath::Pi() + trackAlpha;
561     Double_t secAlpha = 2 * TMath::Pi() / 18.;
562    
563     Int_t sector = static_cast<Int_t>(trackAlpha/secAlpha);
564     Double_t etaTrack = track->Eta();
565     Int_t stack = -1;
566     for(Int_t istack = 0; istack < 5; istack++){
567       if(etaTrack >= etamin[istack] && etaTrack <= etamax[istack]){
568         stack = istack;
569         break;
570       }
571     }
572     if(stack < 0) {
573       AliDebug(2, "Dead Area");
574       continue;
575     }
576
577     detectors[ily] = sector * kNStacks * kNPlanes + stack * kNPlanes + ily;
578     nDet++;
579   }
580   return nDet ? kTRUE : kFALSE;
581 }
582