PHOS/AliPHOSClusterizerv1.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 /* $Id$ */
  17
  18 /* History of cvs commits:
  19  *
  20  * $Log$
  21  * Revision 1.81  2005/09/02 14:32:07  kharlov
  22  * Calibration of raw data
  23  *
  24  * Revision 1.80  2005/08/24 15:31:36  kharlov
  25  * Setting raw digits flag
  26  *
  27  * Revision 1.79  2005/07/25 15:53:53  kharlov
  28  * Read raw data
  29  *
  30  * Revision 1.78  2005/05/28 14:19:04  schutz
  31  * Compilation warnings fixed by T.P.
  32  *
  33  */
  34
  35 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  36 //////////////////////////////////////////////////////////////////////////////
  37 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and
  38 //  unfolds the clusters having several local maxima.
  39 //  Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
  40 //  PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
  41 //  parameters including input digits branch title, thresholds etc.)
  42 //  This TTask is normally called from Reconstructioner, but can as well be used in
  43 //  standalone mode.
  44 // Use Case:
  45 //  root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")
  46 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  47 //               // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
  48 //               // and saves recpoints in branch named "recpointsname" (default = "digitsname")
  49 //  root [1] cl->ExecuteTask()
  50 //               //finds RecPoints in all events stored in galice.root
  51 //  root [2] cl->SetDigitsBranch("digits2")
  52 //               //sets another title for Digitis (input) branch
  53 //  root [3] cl->SetRecPointsBranch("recp2")
  54 //               //sets another title four output branches
  55 //  root [4] cl->SetEmcLocalMaxCut(0.03)
  56 //               //set clusterization parameters
  57 //  root [5] cl->ExecuteTask("deb all time")
  58 //               //once more finds RecPoints options are
  59 //               // deb - print number of found rec points
  60 //               // deb all - print number of found RecPoints and some their characteristics
  61 //               // time - print benchmarking results
  62
  63 // --- ROOT system ---
  64
  65 #include "TMath.h"
  66 #include "TMinuit.h"
  67 #include "TTree.h"
  68 #include "TBenchmark.h"
  69
  70 // --- Standard library ---
  71
  72 // --- AliRoot header files ---
  73 #include "AliLog.h"
  74 #include "AliPHOSGetter.h"
  75 #include "AliPHOSGeometry.h"
  76 #include "AliPHOSClusterizerv1.h"
  77 #include "AliPHOSEmcRecPoint.h"
  78 #include "AliPHOSCpvRecPoint.h"
  79 #include "AliPHOSDigit.h"
  80 #include "AliPHOSDigitizer.h"
  81 #include "AliPHOSCalibrationDB.h"
  82
  83 ClassImp(AliPHOSClusterizerv1)
  84
  85 //____________________________________________________________________________
  86   AliPHOSClusterizerv1::AliPHOSClusterizerv1() : AliPHOSClusterizer()
  87 {
  88   // default ctor (to be used mainly by Streamer)
  89
  90   InitParameters() ;
  91   fDefaultInit = kTRUE ;
  92 }
  93
  94 //____________________________________________________________________________
  95 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName)
  96 :AliPHOSClusterizer(alirunFileName, eventFolderName)
  97 {
  98   // ctor with the indication of the file where header Tree and digits Tree are stored
  99
 100   InitParameters() ;
 101   Init() ;
 102   fDefaultInit = kFALSE ;
 103 }
 104
 105 //____________________________________________________________________________
 106   AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 107 {
 108   // dtor
 109
 110 }
 111 //____________________________________________________________________________
 112 const TString AliPHOSClusterizerv1::BranchName() const
 113 {
 114   return GetName();
 115 }
 116
 117 //____________________________________________________________________________
 118 Float_t  AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId)
 119 {
 120   // Convert digitized amplitude into energy.
 121   // Calibration parameters are taken from calibration data base for raw data,
 122   // or from digitizer parameters for simulated data.
 123
 124   if(fCalibData){
 125     Int_t relId[4];
 126     AliPHOSGetter *gime = AliPHOSGetter::Instance();
 127     gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
 128     Int_t   module = relId[0];
 129     Int_t   column = relId[2];
 130     Int_t   row    = relId[3];
 131     if(absId <= fEmcCrystals) { //calibrate as EMC
 132       fADCchanelEmc   = fCalibData->GetADCchannelEmc (module,column,row);
 133       fADCpedestalEmc = fCalibData->GetADCpedestalEmc(module,column,row);
 134       return fADCpedestalEmc + amp*fADCchanelEmc ;
 135     }
 136     else //calibrate as CPV, not implemented yet
 137       return 0;
 138   }
 139   else{ //simulation
 140     if(absId <= fEmcCrystals) //calibrate as EMC
 141       return fADCpedestalEmc + amp*fADCchanelEmc ;
 142     else //calibrate as CPV
 143       return fADCpedestalCpv+ amp*fADCchanelCpv ;
 144   }
 145 }
 146
 147 //____________________________________________________________________________
 148 void AliPHOSClusterizerv1::Exec(Option_t *option)
 149 {
 150   // Steering method to perform clusterization for events
 151   // in the range from fFirstEvent to fLastEvent.
 152   // This range is optionally set by SetEventRange().
 153   // if fLastEvent=-1 (by default), then process events until the end.
 154
 155   if(strstr(option,"tim"))
 156     gBenchmark->Start("PHOSClusterizer");
 157
 158   if(strstr(option,"print")) {
 159     Print() ;
 160     return ;
 161   }
 162
 163   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 164   if (fRawReader == 0)
 165     gime->SetRawDigits(kFALSE);
 166   else
 167     gime->SetRawDigits(kTRUE);
 168   GetCalibrationParameters() ;
 169
 170   if (fLastEvent == -1)
 171     fLastEvent = gime->MaxEvent() - 1 ;
 172   else
 173     fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // one event at the time
 174   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 175
 176   Int_t ievent ;
 177   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
 178     if (fRawReader == 0)
 179       gime->Event(ievent    ,"D"); // Read digits from simulated data
 180     else
 181       gime->Event(fRawReader,"W"); // Read digits from raw data
 182
 183     fNumberOfEmcClusters  = fNumberOfCpvClusters  = 0 ;
 184
 185     MakeClusters() ;
 186
 187     if(fToUnfold)
 188       MakeUnfolding() ;
 189
 190     WriteRecPoints();
 191
 192     if(strstr(option,"deb"))
 193       PrintRecPoints(option) ;
 194
 195     //increment the total number of recpoints per run
 196     fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
 197     fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
 198     }
 199
 200   if(fWrite) //do not unload in "on flight" mode
 201     Unload();
 202
 203   if(strstr(option,"tim")){
 204     gBenchmark->Stop("PHOSClusterizer");
 205     AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n",
 206          gBenchmark->GetCpuTime("PHOSClusterizer"),
 207          gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ;
 208   }
 209 }
 210
 211 //____________________________________________________________________________
 212 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
 213                                     Int_t nPar, Float_t * fitparameters) const
 214 {
 215   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 216   // The initial values for fitting procedure are set equal to the positions of local maxima.
 217   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 218
 219
 220   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 221   TClonesArray * digits = gime->Digits();
 222
 223
 224   gMinuit->mncler();                     // Reset Minuit's list of paramters
 225   gMinuit->SetPrintLevel(-1) ;           // No Printout
 226   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
 227                                          // To set the address of the minimization function
 228
 229   TList * toMinuit = new TList();
 230   toMinuit->AddAt(emcRP,0) ;
 231   toMinuit->AddAt(digits,1) ;
 232
 233   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 234
 235   // filling initial values for fit parameters
 236   AliPHOSDigit * digit ;
 237
 238   Int_t ierflg  = 0;
 239   Int_t index   = 0 ;
 240   Int_t nDigits = (Int_t) nPar / 3 ;
 241
 242   Int_t iDigit ;
 243
 244   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 245
 246   for(iDigit = 0; iDigit < nDigits; iDigit++){
 247     digit = maxAt[iDigit];
 248
 249     Int_t relid[4] ;
 250     Float_t x = 0.;
 251     Float_t z = 0.;
 252     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 253     geom->RelPosInModule(relid, x, z) ;
 254
 255     Float_t energy = maxAtEnergy[iDigit] ;
 256
 257     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 258     index++ ;
 259     if(ierflg != 0){
 260       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
 261       return kFALSE;
 262     }
 263     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 264     index++ ;
 265     if(ierflg != 0){
 266        Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
 267       return kFALSE;
 268     }
 269     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 270     index++ ;
 271     if(ierflg != 0){
 272       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
 273       return kFALSE;
 274     }
 275   }
 276
 277   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 278                       //  depends on it.
 279   Double_t p1 = 1.0 ;
 280   Double_t p2 = 0.0 ;
 281
 282   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 283   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 284   gMinuit->SetMaxIterations(5);
 285   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 286
 287   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 288
 289   if(ierflg == 4){  // Minimum not found
 290     Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
 291     return kFALSE ;
 292   }
 293   for(index = 0; index < nPar; index++){
 294     Double_t err ;
 295     Double_t val ;
 296     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 297     fitparameters[index] = val ;
 298    }
 299
 300   delete toMinuit ;
 301   return kTRUE;
 302
 303 }
 304
 305 //____________________________________________________________________________
 306 void AliPHOSClusterizerv1::GetCalibrationParameters()
 307 {
 308   // Set calibration parameters:
 309   // For raw data they are read from the calibration database,
 310   // for simulated data they are taken from digitizer.
 311   //
 312   // It is a user responsilibity to open CDB and set
 313   // AliPHOSCalibData object by the following operators:
 314   //
 315   // AliCDBLocal *loc = new AliCDBLocal("deCalibDB");
 316   // AliPHOSCalibData* clb = (AliPHOSCalibData*)AliCDBStorage::Instance()
 317   //    ->Get(path_to_calibdata,run_number);
 318   // AliPHOSGetter* gime = AliPHOSGetter::Instance("galice.root");
 319   // gime->SetCalibData(clb);
 320
 321   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 322   if(gime->IsRawDigits()){
 323     fCalibData = gime->CalibData();
 324   }
 325   else{
 326     if ( !gime->Digitizer() )
 327       gime->LoadDigitizer();
 328     AliPHOSDigitizer * dig = gime->Digitizer();
 329     fADCchanelEmc   = dig->GetEMCchannel() ;
 330     fADCpedestalEmc = dig->GetEMCpedestal();
 331
 332     fADCchanelCpv   = dig->GetCPVchannel() ;
 333     fADCpedestalCpv = dig->GetCPVpedestal() ;
 334   }
 335 }
 336
 337 //____________________________________________________________________________
 338 void AliPHOSClusterizerv1::Init()
 339 {
 340   // Make all memory allocations which can not be done in default constructor.
 341   // Attach the Clusterizer task to the list of PHOS tasks
 342
 343   AliPHOSGetter* gime = AliPHOSGetter::Instance() ;
 344   if(!gime)
 345     gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
 346
 347   AliPHOSGeometry * geom = gime->PHOSGeometry();
 348
 349   fEmcCrystals = geom->GetNModules() *  geom->GetNCristalsInModule() ;
 350
 351   if(!gMinuit)
 352     gMinuit = new TMinuit(100);
 353
 354   if ( !gime->Clusterizer() ) {
 355     gime->PostClusterizer(this);
 356   }
 357 }
 358
 359 //____________________________________________________________________________
 360 void AliPHOSClusterizerv1::InitParameters()
 361 {
 362
 363   fNumberOfCpvClusters     = 0 ;
 364   fNumberOfEmcClusters     = 0 ;
 365
 366   fCpvClusteringThreshold  = 0.0;
 367   fEmcClusteringThreshold  = 0.2;
 368
 369   fEmcLocMaxCut            = 0.03 ;
 370   fCpvLocMaxCut            = 0.03 ;
 371
 372   fEmcMinE                 = 0.01 ;
 373   fCpvMinE                 = 0.0  ;
 374
 375   fW0                      = 4.5 ;
 376   fW0CPV                   = 4.0 ;
 377
 378   fEmcTimeGate             = 1.e-8 ;
 379
 380   fToUnfold                = kTRUE ;
 381
 382   fRecPointsInRun          = 0 ;
 383
 384   fWrite                   = kTRUE ;
 385
 386   fCalibData               = 0 ;
 387
 388   SetEventRange(0,-1) ;
 389 }
 390
 391 //____________________________________________________________________________
 392 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
 393 {
 394   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 395   //                                       = 1 are neighbour
 396   //                                       = 2 are not neighbour but do not continue searching
 397   // neighbours are defined as digits having at least a common vertex
 398   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 399   //                                      which is compared to a digit (d2)  not yet in a cluster
 400
 401   AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
 402
 403   Int_t rv = 0 ;
 404
 405   Int_t relid1[4] ;
 406   geom->AbsToRelNumbering(d1->GetId(), relid1) ;
 407
 408   Int_t relid2[4] ;
 409   geom->AbsToRelNumbering(d2->GetId(), relid2) ;
 410
 411   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
 412     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 413     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 414
 415     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
 416       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
 417       rv = 1 ;
 418     }
 419     else {
 420       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 421         rv = 2; //  Difference in row numbers is too large to look further
 422     }
 423
 424   }
 425   else {
 426
 427     if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
 428       rv=2 ;
 429
 430   }
 431
 432   return rv ;
 433 }
 434 //____________________________________________________________________________
 435 void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits){
 436   for(Int_t i=0; i<digits->GetEntriesFast(); i++){
 437     AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
 438     Float_t cut = IsInEmc(digit) ? fEmcMinE : fCpvMinE ;
 439     if(Calibrate(digit->GetAmp(),digit->GetId()) < cut)
 440       digits->RemoveAt(i) ;
 441   }
 442   digits->Compress() ;
 443   for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
 444     AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 445     digit->SetIndexInList(i) ;
 446   }
 447
 448 }
 449 //____________________________________________________________________________
 450 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 451 {
 452   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
 453
 454   Bool_t rv = kFALSE ;
 455   AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
 456
 457   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 458
 459   if(digit->GetId() <= nEMC )   rv = kTRUE;
 460
 461   return rv ;
 462 }
 463
 464 //____________________________________________________________________________
 465 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 466 {
 467   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
 468
 469   Bool_t rv = kFALSE ;
 470
 471   AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
 472
 473   Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
 474
 475   if(digit->GetId() > nEMC )   rv = kTRUE;
 476
 477   return rv ;
 478 }
 479
 480 //____________________________________________________________________________
 481 void AliPHOSClusterizerv1::Unload()
 482 {
 483   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 484   gime->PhosLoader()->UnloadDigits() ;
 485   gime->PhosLoader()->UnloadRecPoints() ;
 486 }
 487
 488 //____________________________________________________________________________
 489 void AliPHOSClusterizerv1::WriteRecPoints()
 490 {
 491
 492   // Creates new branches with given title
 493   // fills and writes into TreeR.
 494
 495   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 496
 497   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 498   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 499   TClonesArray * digits = gime->Digits() ;
 500
 501
 502   Int_t index ;
 503   //Evaluate position, dispersion and other RecPoint properties..
 504   Int_t nEmc = emcRecPoints->GetEntriesFast();
 505   for(index = 0; index < nEmc; index++){
 506     AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
 507     rp->Purify(fEmcMinE) ;
 508     if(rp->GetMultiplicity()>0.) //If this RP is not empty
 509       rp->EvalAll(fW0,digits) ;
 510     else{
 511       emcRecPoints->RemoveAt(index) ;
 512       delete rp ;
 513     }
 514   }
 515   emcRecPoints->Compress() ;
 516   emcRecPoints->Sort() ;
 517   //  emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
 518   for(index = 0; index < emcRecPoints->GetEntries(); index++){
 519     dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
 520   }
 521
 522   //Now the same for CPV
 523   for(index = 0; index < cpvRecPoints->GetEntries(); index++){
 524     AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
 525     rp->EvalAll(fW0CPV,digits) ;
 526   }
 527   cpvRecPoints->Sort() ;
 528
 529   for(index = 0; index < cpvRecPoints->GetEntries(); index++)
 530     dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
 531
 532   cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
 533
 534   if(fWrite){ //We write TreeR
 535     TTree * treeR = gime->TreeR();
 536
 537     Int_t bufferSize = 32000 ;
 538     Int_t splitlevel = 0 ;
 539
 540     //First EMC
 541     TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
 542     emcBranch->SetTitle(BranchName());
 543
 544     //Now CPV branch
 545     TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
 546     cpvBranch->SetTitle(BranchName());
 547
 548     emcBranch        ->Fill() ;
 549     cpvBranch        ->Fill() ;
 550
 551     gime->WriteRecPoints("OVERWRITE");
 552     gime->WriteClusterizer("OVERWRITE");
 553   }
 554 }
 555
 556 //____________________________________________________________________________
 557 void AliPHOSClusterizerv1::MakeClusters()
 558 {
 559   // Steering method to construct the clusters stored in a list of Reconstructed Points
 560   // A cluster is defined as a list of neighbour digits
 561
 562
 563   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 564
 565   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 566   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 567
 568   emcRecPoints->Delete() ;
 569   cpvRecPoints->Delete() ;
 570
 571   TClonesArray * digits = gime->Digits() ;
 572
 573   //Remove digits below threshold
 574   CleanDigits(digits) ;
 575
 576
 577   TClonesArray * digitsC =  static_cast<TClonesArray*>( digits->Clone() ) ;
 578
 579
 580   // Clusterization starts
 581
 582   TIter nextdigit(digitsC) ;
 583   AliPHOSDigit * digit ;
 584   Bool_t notremoved = kTRUE ;
 585
 586   while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
 587
 588
 589     AliPHOSRecPoint * clu = 0 ;
 590
 591     TArrayI clusterdigitslist(1500) ;
 592     Int_t index ;
 593
 594     if (( IsInEmc (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fEmcClusteringThreshold  ) ||
 595         ( IsInCpv (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fCpvClusteringThreshold  ) ) {
 596       Int_t iDigitInCluster = 0 ;
 597
 598       if  ( IsInEmc(digit) ) {
 599         // start a new EMC RecPoint
 600         if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
 601           emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
 602
 603         emcRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
 604         clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
 605           fNumberOfEmcClusters++ ;
 606         clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId())) ;
 607         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
 608         iDigitInCluster++ ;
 609         digitsC->Remove(digit) ;
 610
 611       } else {
 612
 613         // start a new CPV cluster
 614         if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
 615           cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
 616
 617         cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
 618
 619         clu =  dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
 620         fNumberOfCpvClusters++ ;
 621         clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId()) ) ;
 622         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;
 623         iDigitInCluster++ ;
 624         digitsC->Remove(digit) ;
 625         nextdigit.Reset() ;
 626
 627         // Here we remove remaining EMC digits, which cannot make a cluster
 628
 629         if( notremoved ) {
 630           while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
 631             if( IsInEmc(digit) )
 632               digitsC->Remove(digit) ;
 633             else
 634               break ;
 635           }
 636           notremoved = kFALSE ;
 637         }
 638
 639       } // else
 640
 641       nextdigit.Reset() ;
 642
 643       AliPHOSDigit * digitN ;
 644       index = 0 ;
 645       while (index < iDigitInCluster){ // scan over digits already in cluster
 646         digit =  dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) )  ;
 647         index++ ;
 648         while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
 649           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 650           switch (ineb ) {
 651           case 0 :   // not a neighbour
 652             break ;
 653           case 1 :   // are neighbours
 654             clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->GetId() ) ) ;
 655             clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
 656             iDigitInCluster++ ;
 657             digitsC->Remove(digitN) ;
 658             break ;
 659           case 2 :   // too far from each other
 660             goto endofloop;
 661           } // switch
 662
 663         } // while digitN
 664
 665       endofloop: ;
 666         nextdigit.Reset() ;
 667
 668       } // loop over cluster
 669
 670     } // energy theshold
 671
 672
 673   } // while digit
 674
 675   delete digitsC ;
 676
 677 }
 678
 679 //____________________________________________________________________________
 680 void AliPHOSClusterizerv1::MakeUnfolding()
 681 {
 682   // Unfolds clusters using the shape of an ElectroMagnetic shower
 683   // Performs unfolding of all EMC/CPV clusters
 684
 685   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 686
 687   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 688
 689   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 690   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 691   TClonesArray * digits = gime->Digits() ;
 692
 693   // Unfold first EMC clusters
 694   if(fNumberOfEmcClusters > 0){
 695
 696     Int_t nModulesToUnfold = geom->GetNModules() ;
 697
 698     Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
 699     Int_t index ;
 700     for(index = 0 ; index < numberofNotUnfolded ; index++){
 701
 702       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
 703       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
 704         break ;
 705
 706       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 707       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 708       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 709       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
 710
 711       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 712         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 713         emcRecPoints->Remove(emcRecPoint);
 714         emcRecPoints->Compress() ;
 715         index-- ;
 716         fNumberOfEmcClusters -- ;
 717         numberofNotUnfolded-- ;
 718       }
 719       else{
 720         emcRecPoint->SetNExMax(1) ; //Only one local maximum
 721       }
 722
 723       delete[] maxAt ;
 724       delete[] maxAtEnergy ;
 725     }
 726   }
 727   // Unfolding of EMC clusters finished
 728
 729
 730   // Unfold now CPV clusters
 731   if(fNumberOfCpvClusters > 0){
 732
 733     Int_t nModulesToUnfold = geom->GetNModules() ;
 734
 735     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
 736     Int_t index ;
 737     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
 738
 739       AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
 740
 741       if(recPoint->GetPHOSMod()> nModulesToUnfold)
 742         break ;
 743
 744       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
 745
 746       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 747       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 748       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 749       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
 750
 751       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 752         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 753         cpvRecPoints->Remove(emcRecPoint);
 754         cpvRecPoints->Compress() ;
 755         index-- ;
 756         numberofCpvNotUnfolded-- ;
 757         fNumberOfCpvClusters-- ;
 758       }
 759
 760       delete[] maxAt ;
 761       delete[] maxAtEnergy ;
 762     }
 763   }
 764   //Unfolding of Cpv clusters finished
 765
 766 }
 767
 768 //____________________________________________________________________________
 769 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t r)
 770 {
 771   // Shape of the shower (see PHOS TDR)
 772   // If you change this function, change also the gradient evaluation in ChiSquare()
 773
 774   Double_t r4    = r*r*r*r ;
 775   Double_t r295  = TMath::Power(r, 2.95) ;
 776   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 777   return shape ;
 778 }
 779
 780 //____________________________________________________________________________
 781 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 782                                           Int_t nMax,
 783                                           AliPHOSDigit ** maxAt,
 784                                           Float_t * maxAtEnergy)
 785 {
 786   // Performs the unfolding of a cluster with nMax overlapping showers
 787
 788   AliPHOSGetter * gime = AliPHOSGetter::Instance();
 789
 790   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 791
 792   const TClonesArray * digits = gime->Digits() ;
 793   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
 794   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
 795
 796   Int_t nPar = 3 * nMax ;
 797   Float_t * fitparameters = new Float_t[nPar] ;
 798
 799   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 800   if( !rv ) {
 801     // Fit failed, return and remove cluster
 802     iniEmc->SetNExMax(-1) ;
 803     delete[] fitparameters ;
 804     return ;
 805   }
 806
 807   // create ufolded rec points and fill them with new energy lists
 808   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
 809   // and later correct this number in acordance with actual energy deposition
 810
 811   Int_t nDigits = iniEmc->GetMultiplicity() ;
 812   Float_t * efit = new Float_t[nDigits] ;
 813   Float_t xDigit=0.,zDigit=0.,distance=0. ;
 814   Float_t xpar=0.,zpar=0.,epar=0.  ;
 815   Int_t relid[4] ;
 816   AliPHOSDigit * digit = 0 ;
 817   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 818
 819   Int_t iparam ;
 820   Int_t iDigit ;
 821   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 822     digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;
 823     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 824     geom->RelPosInModule(relid, xDigit, zDigit) ;
 825     efit[iDigit] = 0;
 826
 827     iparam = 0 ;
 828     while(iparam < nPar ){
 829       xpar = fitparameters[iparam] ;
 830       zpar = fitparameters[iparam+1] ;
 831       epar = fitparameters[iparam+2] ;
 832       iparam += 3 ;
 833       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 834       distance =  TMath::Sqrt(distance) ;
 835       efit[iDigit] += epar * ShowerShape(distance) ;
 836     }
 837   }
 838
 839
 840   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
 841   // so that energy deposited in each cell is distributed betwin new clusters proportionally
 842   // to its contribution to efit
 843
 844   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 845   Float_t ratio ;
 846
 847   iparam = 0 ;
 848   while(iparam < nPar ){
 849     xpar = fitparameters[iparam] ;
 850     zpar = fitparameters[iparam+1] ;
 851     epar = fitparameters[iparam+2] ;
 852     iparam += 3 ;
 853
 854     AliPHOSEmcRecPoint * emcRP = 0 ;
 855
 856     if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
 857
 858       if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
 859         emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
 860
 861       (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
 862       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
 863       fNumberOfEmcClusters++ ;
 864       emcRP->SetNExMax((Int_t)nPar/3) ;
 865     }
 866     else{//create new entries in fCpvRecPoints
 867       if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
 868         cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
 869
 870       (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
 871       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
 872       fNumberOfCpvClusters++ ;
 873     }
 874
 875     Float_t eDigit ;
 876     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 877       digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ;
 878       geom->AbsToRelNumbering(digit->GetId(), relid) ;
 879       geom->RelPosInModule(relid, xDigit, zDigit) ;
 880       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 881       distance =  TMath::Sqrt(distance) ;
 882       ratio = epar * ShowerShape(distance) / efit[iDigit] ;
 883       eDigit = emcEnergies[iDigit] * ratio ;
 884       emcRP->AddDigit( *digit, eDigit ) ;
 885     }
 886   }
 887
 888   delete[] fitparameters ;
 889   delete[] efit ;
 890
 891 }
 892
 893 //_____________________________________________________________________________
 894 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 895 {
 896   // Calculates the Chi square for the cluster unfolding minimization
 897   // Number of parameters, Gradient, Chi squared, parameters, what to do
 898
 899   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
 900
 901   AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
 902   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
 903
 904
 905
 906   //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
 907
 908   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 909
 910   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
 911
 912   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
 913
 914   const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
 915   fret = 0. ;
 916   Int_t iparam ;
 917
 918   if(iflag == 2)
 919     for(iparam = 0 ; iparam < nPar ; iparam++)
 920       Grad[iparam] = 0 ; // Will evaluate gradient
 921
 922   Double_t efit ;
 923
 924   AliPHOSDigit * digit ;
 925   Int_t iDigit ;
 926
 927   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 928
 929     digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
 930
 931     Int_t relid[4] ;
 932     Float_t xDigit ;
 933     Float_t zDigit ;
 934
 935     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 936
 937     geom->RelPosInModule(relid, xDigit, zDigit) ;
 938
 939      if(iflag == 2){  // calculate gradient
 940        Int_t iParam = 0 ;
 941        efit = 0 ;
 942        while(iParam < nPar ){
 943          Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
 944          iParam++ ;
 945          distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
 946          distance = TMath::Sqrt( distance ) ;
 947          iParam++ ;
 948          efit += x[iParam] * ShowerShape(distance) ;
 949          iParam++ ;
 950        }
 951        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
 952        iParam = 0 ;
 953        while(iParam < nPar ){
 954          Double_t xpar = x[iParam] ;
 955          Double_t zpar = x[iParam+1] ;
 956          Double_t epar = x[iParam+2] ;
 957          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
 958          Double_t shape = sum * ShowerShape(dr) ;
 959          Double_t r4 = dr*dr*dr*dr ;
 960          Double_t r295 = TMath::Power(dr,2.95) ;
 961          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
 962                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
 963
 964          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
 965          iParam++ ;
 966          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
 967          iParam++ ;
 968          Grad[iParam] += shape ;                                  // Derivative over energy
 969          iParam++ ;
 970        }
 971      }
 972      efit = 0;
 973      iparam = 0 ;
 974
 975      while(iparam < nPar ){
 976        Double_t xpar = x[iparam] ;
 977        Double_t zpar = x[iparam+1] ;
 978        Double_t epar = x[iparam+2] ;
 979        iparam += 3 ;
 980        Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
 981        distance =  TMath::Sqrt(distance) ;
 982        efit += epar * ShowerShape(distance) ;
 983      }
 984
 985      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
 986      // Here we assume, that sigma = sqrt(E)
 987   }
 988
 989 }
 990
 991 //____________________________________________________________________________
 992 void AliPHOSClusterizerv1::Print(const Option_t *)const
 993 {
 994   // Print clusterizer parameters
 995
 996   TString message ;
 997   TString taskName(GetName()) ;
 998   taskName.ReplaceAll(Version(), "") ;
 999
1000   if( strcmp(GetName(), "") !=0 ) {
1001     // Print parameters
1002     message  = "\n--------------- %s %s -----------\n" ;
1003     message += "Clusterizing digits from the file: %s\n" ;
1004     message += "                           Branch: %s\n" ;
1005     message += "                       EMC Clustering threshold = %f\n" ;
1006     message += "                       EMC Local Maximum cut    = %f\n" ;
1007     message += "                       EMC Logarothmic weight   = %f\n" ;
1008     message += "                       CPV Clustering threshold = %f\n" ;
1009     message += "                       CPV Local Maximum cut    = %f\n" ;
1010     message += "                       CPV Logarothmic weight   = %f\n" ;
1011     if(fToUnfold)
1012       message += " Unfolding on\n" ;
1013     else
1014       message += " Unfolding off\n" ;
1015
1016     message += "------------------------------------------------------------------" ;
1017   }
1018   else
1019     message = " AliPHOSClusterizerv1 not initialized " ;
1020
1021   AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
1022        taskName.Data(),
1023        GetTitle(),
1024        taskName.Data(),
1025        GetName(),
1026        fEmcClusteringThreshold,
1027        fEmcLocMaxCut,
1028        fW0,
1029        fCpvClusteringThreshold,
1030        fCpvLocMaxCut,
1031        fW0CPV )) ;
1032 }
1033
1034
1035 //____________________________________________________________________________
1036 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1037 {
1038   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1039
1040   AliPHOSGetter * gime = AliPHOSGetter::Instance();
1041
1042   TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1043   TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
1044
1045   AliInfo(Form("\nevent %d \n    Found %d EMC RecPoints and %d CPV RecPoints",
1046                gAlice->GetEvNumber(),
1047                emcRecPoints->GetEntriesFast(),
1048                cpvRecPoints->GetEntriesFast() ))  ;
1049
1050   fRecPointsInRun +=  emcRecPoints->GetEntriesFast() ;
1051   fRecPointsInRun +=  cpvRecPoints->GetEntriesFast() ;
1052
1053
1054   if(strstr(option,"all")) {
1055     printf("\n EMC clusters \n") ;
1056     printf("Index    Ene(MeV) Multi Module     X    Y   Z    Lambdas_1  Lambda_2  # of prim  Primaries list\n") ;
1057     Int_t index ;
1058     for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
1059       AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
1060       TVector3  locpos;
1061       rp->GetLocalPosition(locpos);
1062       Float_t lambda[2];
1063       rp->GetElipsAxis(lambda);
1064       Int_t * primaries;
1065       Int_t nprimaries;
1066       primaries = rp->GetPrimaries(nprimaries);
1067       printf("\n%6d  %8.2f  %3d     %2d     %4.1f   %4.1f   %4.1f   %4f  %4f    %2d     : ",
1068               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1069               locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1070
1071       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1072         printf("%d ", primaries[iprimary] ) ;
1073       }
1074       printf("\n") ;
1075     }
1076
1077     //Now plot CPV recPoints
1078     printf("\n CPV clusters \n") ;
1079     printf("Index    Ene(MeV) Module     X     Y    Z  \n") ;
1080     for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
1081       AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )cpvRecPoints->At(index) ;
1082
1083       TVector3  locpos;
1084       rp->GetLocalPosition(locpos);
1085
1086       printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f \n",
1087              rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1088              locpos.X(), locpos.Y(), locpos.Z()) ;
1089     }
1090   }
1091 }
1092