FMD/AliFMDBackgroundCorrection.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2004, 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 // Thil class computes background corrections for the FMD. The correction is computed
  17 // in eta,phi cells and the objects stored can be put into alien to use with analysis.
  18 // It is based on the AliFMDInput class that is used to loop over hits and primaries.
  19 //
  20 // Author: Hans Hjersing Dalsgaard, NBI, hans.dalsgaard@cern.ch
  21 //
  22 //
  23
  24 #include "AliSimulation.h"
  25 #include "TStopwatch.h"
  26 #include "iostream"
  27 #include "TGrid.h"
  28 #include "AliRunLoader.h"
  29 #include "AliGeomManager.h"
  30 #include "AliFMDGeometry.h"
  31 #include "AliStack.h"
  32 #include "TParticle.h"
  33 #include "TDatabasePDG.h"
  34 #include "TParticlePDG.h"
  35 #include "AliRun.h"
  36 #include "AliFMDBackgroundCorrection.h"
  37 #include "TSystem.h"
  38 #include "AliCDBManager.h"
  39 #include "TTree.h"
  40 #include "TClonesArray.h"
  41 #include "TBranch.h"
  42 #include "AliFMDHit.h"
  43 #include "AliLoader.h"
  44 #include "AliFMD.h"
  45 #include "AliGenEventHeader.h"
  46 #include "AliHeader.h"
  47 #include "TFile.h"
  48 #include "TAxis.h"
  49 #include "AliCDBId.h"
  50 #include "AliCDBMetaData.h"
  51 #include "TROOT.h"
  52 #include "AliFMDParameters.h"
  53 #include "AliLog.h"
  54 #include "TList.h"
  55 #include "AliFMDAnaParameters.h"
  56 #include "AliFMDAnaCalibBackgroundCorrection.h"
  57 #include "AliTrackReference.h"
  58 #include "AliFMDStripIndex.h"
  59
  60 ClassImp(AliFMDBackgroundCorrection)
  61 //_____________________________________________________________________
  62 AliFMDBackgroundCorrection::AliFMDBackgroundCorrection() :
  63   TNamed(),
  64   fCorrectionArray(),
  65   fPrimaryList()
  66 {}
  67
  68 //_____________________________________________________________________
  69 AliFMDBackgroundCorrection::AliFMDInputBG::AliFMDInputBG(Bool_t hits_not_trackref) :
  70   AliFMDInput("galice.root"),
  71   fPrimaryArray(),
  72   fHitArray(),
  73   fPrimaryMapInner(),
  74   fPrimaryMapOuter(),
  75   fHitMap(),
  76   fLastTrackByStrip(),
  77   fPrim(0),
  78   fHits(0),
  79   fZvtxCut(0),
  80   fNvtxBins(0),
  81   fPrevTrack(-1),
  82   fPrevDetector(-1),
  83   fPrevRing('Q'),
  84   fPrevSec(-1),
  85   fNbinsEta(100)
  86 {
  87
  88   if(hits_not_trackref) {
  89     AddLoad(kHits);
  90     AddLoad(kTracks);
  91   }
  92   else
  93     AddLoad(kTrackRefs);
  94   AddLoad(kKinematics);
  95   AddLoad(kHeader);
  96
  97
  98 }
  99
 100 //_____________________________________________________________________
 101
 102 void AliFMDBackgroundCorrection::GenerateBackgroundCorrection(Bool_t from_hits,
 103                                                               const Int_t nvtxbins,
 104                                                               Float_t zvtxcut,
 105                                                               const Int_t nBinsEta,
 106                                                               Bool_t storeInAlien,
 107                                                               Int_t runNo,
 108                                                               Int_t endRunNo,
 109                                                               const Char_t* filename,
 110                                                               Bool_t simulate,
 111                                                               Int_t nEvents,
 112                                                               Bool_t inFile,
 113                                                               const Char_t* infilename) {
 114
 115   //TGrid::Connect("alien:",0,0,"t");
 116   if(simulate)
 117     Simulate(nEvents);
 118   else {
 119     //AliCDBManager::Instance()->SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
 120     AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 121     AliCDBManager::Instance()->SetRun(runNo);
 122
 123 #if defined(__CINT__)
 124     gSystem->Load("liblhapdf");
 125     gSystem->Load("libEGPythia6");
 126     gSystem->Load("libpythia6");
 127     gSystem->Load("libAliPythia6");
 128     gSystem->Load("libgeant321");
 129 #endif
 130       //
 131   }
 132
 133   //Setting up the geometry
 134   //-----------------------------------------------
 135   if (AliGeomManager::GetGeometry() == NULL)
 136     AliGeomManager::LoadGeometry();
 137
 138   AliFMDGeometry* geo = AliFMDGeometry::Instance();
 139   geo->Init();
 140   geo->InitTransformations();
 141
 142   AliInfo("Processing hits and primaries ");
 143
 144   AliFMDInputBG input(from_hits);
 145
 146   if(!inFile) {
 147
 148     input.SetVtxCutZ(zvtxcut);
 149     input.SetNvtxBins(nvtxbins);
 150     input.SetNbinsEta(nBinsEta);
 151     input.Run();
 152   }
 153
 154   AliInfo(Form("Found %d primaries and %d hits.", input.GetNprim(),input.GetNhits()));
 155   TObjArray* hitArray ;
 156   TObjArray* primaryArray;
 157   if(inFile) {
 158     TFile* infile = TFile::Open(infilename);
 159     hitArray     = new TObjArray();
 160     primaryArray = new TObjArray();
 161
 162     for(Int_t det =1; det<=3;det++)
 163       {
 164         TObjArray* detArrayHits = new TObjArray();
 165         detArrayHits->SetName(Form("FMD%d",det));
 166         hitArray->AddAtAndExpand(detArrayHits,det);
 167         Int_t nRings = (det==1 ? 1 : 2);
 168         for(Int_t ring = 0;ring<nRings;ring++)
 169           {
 170
 171             Char_t ringChar = (ring == 0 ? 'I' : 'O');
 172             TObjArray* vtxArrayHits = new TObjArray();
 173             vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
 174             detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
 175             for(Int_t v=0; v<nvtxbins;v++)
 176               {
 177
 178                 TH2F* hHits          = (TH2F*)infile->Get(Form("hHits_FMD%d%c_vtx%d",det,ringChar,v));
 179
 180
 181               vtxArrayHits->AddAtAndExpand(hHits,v);
 182
 183             }
 184         }
 185
 186     }
 187
 188     for(Int_t iring = 0; iring<2;iring++) {
 189       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 190       TObjArray* ringArray = new TObjArray();
 191       ringArray->SetName(Form("FMD_%c",ringChar));
 192       primaryArray->AddAtAndExpand(ringArray,iring);
 193       for(Int_t v=0; v<nvtxbins;v++) {
 194
 195         TH2F* hPrimary       = (TH2F*)infile->Get(Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
 196         ringArray->AddAtAndExpand(hPrimary,v);
 197       }
 198     }
 199
 200
 201   }
 202   else {
 203     hitArray     = input.GetHits();
 204     primaryArray = input.GetPrimaries();
 205   }
 206   fCorrectionArray.SetName("FMD_bg_correction");
 207   fCorrectionArray.SetOwner();
 208
 209   TList* primaryList     = new TList();
 210   primaryList->SetName("primaries");
 211
 212   TList* hitList     = new TList();
 213   hitList->SetName("hits");
 214   TList* corrList    = new TList();
 215   corrList->SetName("corrections");
 216
 217   AliFMDAnaCalibBackgroundCorrection* background = new AliFMDAnaCalibBackgroundCorrection();
 218
 219   for(Int_t det= 1; det <=3; det++) {
 220     Int_t nRings = (det==1 ? 1 : 2);
 221
 222     TObjArray* detArrayCorrection = new TObjArray();
 223     detArrayCorrection->SetName(Form("FMD%d",det));
 224     fCorrectionArray.AddAtAndExpand(detArrayCorrection,det);
 225
 226
 227     for(Int_t iring = 0; iring<nRings; iring++) {
 228       TObjArray* primRingArray = (TObjArray*)primaryArray->At(iring);
 229       Char_t ringChar = (iring == 0 ? 'I' : 'O');
 230       TObjArray* vtxArrayCorrection = new TObjArray();
 231       vtxArrayCorrection->SetName(Form("FMD%d%c",det,ringChar));
 232       detArrayCorrection->AddAtAndExpand(vtxArrayCorrection,iring);
 233
 234       for(Int_t vertexBin=0;vertexBin<nvtxbins;vertexBin++) {
 235         TObjArray* detArray  = (TObjArray*)hitArray->At(det);
 236         TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
 237         TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
 238         hitList->Add(hHits);
 239         TH2F* hPrimary  = (TH2F*)primRingArray->At(vertexBin);
 240         primaryList->Add(hPrimary);
 241         TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ringChar,vertexBin));
 242         hCorrection->Divide(hPrimary);
 243         corrList->Add(hCorrection);
 244         hCorrection->SetTitle(hCorrection->GetName());
 245         vtxArrayCorrection->AddAtAndExpand(hCorrection,vertexBin);
 246         background->SetBgCorrection(det,ringChar,vertexBin,hCorrection);
 247       }
 248
 249     }
 250   }
 251
 252   TAxis refAxis(nvtxbins,-1*zvtxcut,zvtxcut);
 253   if(inFile) {
 254     TFile* infile = TFile::Open(infilename);
 255     TAxis* refaxis = (TAxis*)infile->Get("vertexbins");
 256     background->SetRefAxis(refaxis);
 257
 258   }
 259   else background->SetRefAxis(&refAxis);
 260
 261   TFile*  fout = new TFile(filename,"RECREATE");
 262   refAxis.Write("vertexbins");
 263
 264   hitList->Write();
 265   primaryList->Write();
 266   corrList->Write();
 267
 268   TObjArray* container = new TObjArray();
 269   container->SetOwner();
 270   container->AddAtAndExpand(&refAxis,0);
 271   container->AddAtAndExpand(&fCorrectionArray,1);
 272   container->AddAtAndExpand(hitArray,2);
 273   container->AddAtAndExpand(primaryArray,3);
 274
 275
 276   if(storeInAlien) {
 277     AliCDBManager* cdb = AliCDBManager::Instance();
 278     cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
 279     AliCDBId      id(AliFMDAnaParameters::GetBackgroundPath(),runNo,endRunNo);
 280
 281     AliCDBMetaData* meta = new AliCDBMetaData;
 282     meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
 283     meta->SetAliRootVersion(gROOT->GetVersion());
 284     meta->SetBeamPeriod(1);
 285     meta->SetComment("Background Correction for FMD");
 286     meta->SetProperty("key1", background );
 287     cdb->Put(background, id, meta);
 288
 289   }
 290
 291   fout->Close();
 292
 293
 294   }
 295
 296 //_____________________________________________________________________
 297 void AliFMDBackgroundCorrection::Simulate(Int_t nEvents) {
 298
 299   AliSimulation sim ;
 300   sim.SetRunNumber(0);
 301   TGrid::Connect("alien:",0,0,"t");
 302   sim.SetDefaultStorage("alien://Folder=/alice/data/2008/LHC08d/OCDB/");
 303   sim.SetConfigFile("Config.C");
 304   sim.SetRunQA("FMD:");
 305   TStopwatch timer;
 306   timer.Start();
 307   sim.RunSimulation(nEvents);
 308   timer.Stop();
 309   timer.Print();
 310
 311 }
 312
 313 //_____________________________________________________________________
 314 Bool_t
 315 AliFMDBackgroundCorrection::AliFMDInputBG::ProcessHit(AliFMDHit* h,
 316                                                       TParticle* /*p*/)
 317 {
 318
 319   if(!h)
 320     return kTRUE;
 321   Bool_t retval = ProcessEvent(h->Detector(),
 322                                h->Ring(),
 323                                h->Sector(),
 324                                h->Strip(),
 325                                h->Track(),
 326                                h->Q());
 327
 328   // std::cout<<"Length: "<<h->Length()<<std::endl;
 329   // std::cout<<"Is stopped "<<h->IsStop()<<"    "<<kTRUE<<std::endl;
 330   return retval;
 331 }
 332 //_____________________________________________________________________
 333 Bool_t
 334 AliFMDBackgroundCorrection::AliFMDInputBG::ProcessTrackRef(AliTrackReference* tr, TParticle* p)
 335 {
 336   if(!tr)
 337     return kTRUE;
 338   UShort_t det,sec,strip;
 339   Char_t   ring;
 340   AliFMDStripIndex::Unpack(tr->UserId(),det,ring,sec,strip);
 341   Int_t         nTrack  = tr->GetTrack();
 342   TDatabasePDG* pdgDB   = TDatabasePDG::Instance();
 343   TParticlePDG* pdgPart = pdgDB->GetParticle(p->GetPdgCode());
 344   Float_t       charge  = (pdgPart ? pdgPart->Charge() : 0);
 345   Bool_t        retval  = ProcessEvent(det,ring,sec,strip,nTrack,charge);
 346   return retval;
 347
 348 }
 349 //_____________________________________________________________________
 350 Bool_t
 351 AliFMDBackgroundCorrection::AliFMDInputBG::ProcessEvent(UShort_t det,
 352                                                         Char_t   ring,
 353                                                         UShort_t sec,
 354                                                         UShort_t strip,
 355                                                         Int_t    nTrack,
 356                                                         Float_t  charge)
 357 {
 358
 359
 360
 361   /*  if(charge !=  0 &&
 362      ((nTrack != fPrevTrack) ||
 363       (det != fPrevDetector))){ ||
 364       (ring != fPrevRing))){    ||
 365                                   (sec != fPrevSec))) {*/
 366   /*Float_t nstrips = (ring =='O' ? 256 : 512);
 367
 368
 369   Float_t prevStripTrack = -1;
 370   if(strip >0)
 371     prevStripTrack = fLastTrackByStrip.operator()(det,ring,sec,strip-1);
 372
 373   Float_t nextStripTrack = -1;
 374   if(strip < (nstrips - 1))
 375     nextStripTrack = fLastTrackByStrip.operator()(det,ring,sec,strip+1);
 376   */
 377
 378   Float_t thisStripTrack = fLastTrackByStrip.operator()(det,ring,0,0);
 379   // if(nTrack == fLastTrackByStrip.operator()(det,ring,sec,strip))
 380   //  std::cout<<"Track # "<<nTrack<<"  failed the cut in "<<det<<"   "<<ring<<"   "<<sec<<"  "<<strip<<std::endl;
 381   // else
 382   //  std::cout<<"Track # "<<nTrack<<"  passed the cut in "<<det<<"   "<<ring<<"   "<<sec<<"  "<<strip<<std::endl;
 383   if(charge != 0 && nTrack != thisStripTrack) {
 384
 385     fHitMap.operator()(det,ring,sec,strip) += 1;
 386     fHits++;
 387     Float_t nstrips = (ring =='O' ? 256 : 512);
 388     fLastTrackByStrip.operator()(det,ring,sec,strip) = (Float_t)nTrack;
 389     if(strip >0)
 390       fLastTrackByStrip.operator()(det,ring,sec,strip-1) = (Float_t)nTrack;
 391     if(strip < (nstrips - 1))
 392       fLastTrackByStrip.operator()(det,ring,sec,strip+1) = (Float_t)nTrack;
 393
 394     fPrevDetector = det;
 395     fPrevRing     = ring;
 396     fPrevSec      = sec;
 397     fPrevTrack    = nTrack;
 398   }
 399
 400
 401
 402   return kTRUE;
 403
 404 }
 405
 406 //_____________________________________________________________________
 407 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Init()
 408 {
 409   fLastTrackByStrip.Reset(-1);
 410
 411   fPrimaryArray.SetOwner();
 412   fPrimaryArray.SetName("FMD_primary");
 413
 414   fPrimaryMapInner.SetBins(fNbinsEta, -6,6, 20, 0, 2*TMath::Pi());
 415   fPrimaryMapOuter.SetBins(fNbinsEta, -6,6, 40, 0, 2*TMath::Pi());
 416   fPrimaryMapInner.SetName("fPrimaryMapInner");
 417   fPrimaryMapInner.SetName("fPrimaryMapOuter");
 418
 419   fPrimaryMapInner.Sumw2();
 420   fPrimaryMapOuter.Sumw2();
 421   for(Int_t iring = 0; iring<2;iring++) {
 422     Char_t ringChar = (iring == 0 ? 'I' : 'O');
 423     TObjArray* ringArray = new TObjArray();
 424     ringArray->SetName(Form("FMD_%c",ringChar));
 425     fPrimaryArray.AddAtAndExpand(ringArray,iring);
 426     Int_t nSec = (iring == 1 ? 40 : 20);
 427     for(Int_t v=0; v<fNvtxBins;v++) {
 428
 429       TH2F* hPrimary       = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
 430                                       Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
 431                                       fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
 432       hPrimary->Sumw2();
 433       ringArray->AddAtAndExpand(hPrimary,v);
 434     }
 435   }
 436
 437
 438   fHitArray.SetOwner();
 439   fHitArray.SetName("FMD_hits");
 440
 441   for(Int_t det =1; det<=3;det++) {
 442     TObjArray* detArrayHits = new TObjArray();
 443     detArrayHits->SetName(Form("FMD%d",det));
 444     fHitArray.AddAtAndExpand(detArrayHits,det);
 445     Int_t nRings = (det==1 ? 1 : 2);
 446     for(Int_t ring = 0;ring<nRings;ring++) {
 447       Int_t nSec = (ring == 1 ? 40 : 20);
 448       Char_t ringChar = (ring == 0 ? 'I' : 'O');
 449       TObjArray* vtxArrayHits = new TObjArray();
 450       vtxArrayHits->SetName(Form("FMD%d%c",det,ringChar));
 451       detArrayHits->AddAtAndExpand(vtxArrayHits,ring);
 452       for(Int_t v=0; v<fNvtxBins;v++) {
 453         TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
 454                                Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
 455                                fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
 456         hHits->Sumw2();
 457         vtxArrayHits->AddAtAndExpand(hHits,v);
 458
 459       }
 460     }
 461   }
 462
 463   AliFMDInput::Init();
 464
 465   return kTRUE;
 466 }
 467 //_____________________________________________________________________
 468
 469 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::Begin(Int_t event )
 470 {
 471
 472   Bool_t             retVal    = AliFMDInput::Begin(event);
 473   AliStack*          partStack = fLoader->Stack();
 474   Int_t              nTracks   = partStack->GetNtrack();
 475
 476
 477   for(Int_t j=0;j<nTracks;j++) {
 478     TParticle* p           = partStack->Particle(j);
 479     TDatabasePDG* pdgDB    = TDatabasePDG::Instance();
 480     TParticlePDG* pdgPart  = pdgDB->GetParticle(p->GetPdgCode());
 481     Float_t       charge   = (pdgPart ? pdgPart->Charge() : 0);
 482     Float_t       phi      = TMath::ATan2(p->Py(),p->Px());
 483
 484     if(phi<0) phi = phi+2*TMath::Pi();
 485     // if(p->Theta() == 0) continue;
 486     Float_t eta   = p->Eta();
 487
 488     // std::cout<<-1*TMath::Log(TMath::Tan(0.5*p->Theta()))<<std::endl;
 489
 490     Bool_t primary = partStack->IsPhysicalPrimary(j);
 491     //(charge!=0)&&(TMath::Abs(p->Vx() - vertex.At(0))<0.01)&&(TMath::Abs(p->Vy() - vertex.At(1))<0.01)&&(TMath::Abs(p->Vz() - vertex.At(2))<0.01);
 492     if(primary && charge !=0) {
 493
 494       fPrim++;
 495
 496       fPrimaryMapInner.Fill(eta,phi);
 497       fPrimaryMapOuter.Fill(eta,phi);
 498     }
 499   }
 500
 501   return retVal;
 502 }
 503 //_____________________________________________________________________
 504 Bool_t AliFMDBackgroundCorrection::AliFMDInputBG::End()  {
 505
 506   Bool_t retval = AliFMDInput::End();
 507
 508   AliGenEventHeader* genHeader = fLoader->GetHeader()->GenEventHeader();
 509   TArrayF vertex;
 510   genHeader->PrimaryVertex(vertex);
 511
 512   if(TMath::Abs(vertex.At(2)) > fZvtxCut)
 513     return kTRUE;
 514
 515   Double_t delta           = 2*fZvtxCut/fNvtxBins;
 516   Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
 517   Int_t    vertexBin       = (Int_t)vertexBinDouble;
 518   //Primaries
 519   TObjArray* innerArray = (TObjArray*)fPrimaryArray.At(0);
 520   TObjArray* outerArray = (TObjArray*)fPrimaryArray.At(1);
 521
 522   TH2F* hPrimaryInner  = (TH2F*)innerArray->At(vertexBin);
 523   TH2F* hPrimaryOuter  = (TH2F*)outerArray->At(vertexBin);
 524
 525   hPrimaryInner->Add(&fPrimaryMapInner);
 526   hPrimaryOuter->Add(&fPrimaryMapOuter);
 527
 528   //Hits
 529   for(UShort_t det=1;det<=3;det++) {
 530     Int_t nRings = (det==1 ? 1 : 2);
 531     for (UShort_t ir = 0; ir < nRings; ir++) {
 532       Char_t   ring = (ir == 0 ? 'I' : 'O');
 533       UShort_t nsec = (ir == 0 ? 20  : 40);
 534       UShort_t nstr = (ir == 0 ? 512 : 256);
 535
 536       for(UShort_t sec =0; sec < nsec;  sec++) {
 537
 538         for(UShort_t strip = 0; strip < nstr; strip++) {
 539
 540           if(fHitMap.operator()(det,ring,sec,strip) > 0) {
 541             //std::cout<<fHitMap.operator()(det,ring,sec,strip)<<std::endl;
 542             Double_t x,y,z;
 543             AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
 544             fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
 545
 546             Int_t iring = (ring == 'I' ? 0 : 1);
 547
 548             TObjArray* detArray  = (TObjArray*)fHitArray.At(det);
 549             TObjArray* vtxArray  = (TObjArray*)detArray->At(iring);
 550             TH2F* hHits          = (TH2F*)vtxArray->At(vertexBin);
 551
 552             Float_t   phi   = TMath::ATan2(y,x);
 553             if(phi<0) phi   = phi+2*TMath::Pi();
 554             Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
 555             Float_t   theta = TMath::ATan2(r,z-vertex.At(2));
 556             Float_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 557             hHits->Fill(eta,phi,fHitMap.operator()(det,ring,sec,strip));
 558           }
 559
 560         }
 561       }
 562     }
 563   }
 564
 565   fPrimaryMapInner.Reset();
 566   fPrimaryMapOuter.Reset();
 567   fHitMap.Reset(0);
 568   fLastTrackByStrip.Reset(-1);
 569
 570   return retval;
 571 }
 572
 573 //
 574 // EOF
 575 //
 576
 577