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.115  2007/09/26 14:22:17  cvetan
  22  * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
  23  *
  24  * Revision 1.114  2007/09/06 16:06:44  kharlov
  25  * Absence of sorting results in loose of all unfolded clusters
  26  *
  27  * Revision 1.113  2007/08/28 12:55:07  policheh
  28  * Loaders removed from the reconstruction code (C.Cheshkov)
  29  *
  30  * Revision 1.112  2007/08/22 09:20:50  hristov
  31  * Updated QA classes (Yves)
  32  *
  33  * Revision 1.111  2007/08/08 12:11:28  kharlov
  34  * Protection against uninitialized fQADM
  35  *
  36  * Revision 1.110  2007/08/07 14:16:00  kharlov
  37  * Quality assurance added (Yves Schutz)
  38  *
  39  * Revision 1.109  2007/07/24 17:20:35  policheh
  40  * Usage of RecoParam objects instead of hardcoded parameters in reconstruction.
  41  * (See $ALICE_ROOT/PHOS/macros/BeamTest2006/RawReconstruction.C).
  42  *
  43  * Revision 1.108  2007/06/18 07:00:51  kharlov
  44  * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
  45  *
  46  * Revision 1.107  2007/05/25 14:12:26  policheh
  47  * Local to tracking CS transformation added for CPV rec. points
  48  *
  49  * Revision 1.106  2007/05/24 13:01:22  policheh
  50  * Local to tracking CS transformation invoked for each EMC rec.point
  51  *
  52  * Revision 1.105  2007/05/02 13:41:22  kharlov
  53  * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
  54  *
  55  * Revision 1.104  2007/04/27 16:55:53  kharlov
  56  * Calibration stops if PHOS CDB objects do not exist
  57  *
  58  * Revision 1.103  2007/04/11 11:55:45  policheh
  59  * SetDistancesToBadChannels() added.
  60  *
  61  * Revision 1.102  2007/03/28 19:18:15  kharlov
  62  * RecPoints recalculation in TSM removed
  63  *
  64  * Revision 1.101  2007/03/06 06:51:27  kharlov
  65  * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
  66  *
  67  * Revision 1.100  2007/01/10 11:07:26  kharlov
  68  * Raw digits writing to file (B.Polichtchouk)
  69  *
  70  * Revision 1.99  2006/11/07 16:49:51  kharlov
  71  * Corrections for next event switch in case of raw data (B.Polichtchouk)
  72  *
  73  * Revision 1.98  2006/10/27 17:14:27  kharlov
  74  * Introduce AliDebug and AliLog (B.Polichtchouk)
  75  *
  76  * Revision 1.97  2006/08/29 11:41:19  kharlov
  77  * Missing implementation of ctors and = operator are added
  78  *
  79  * Revision 1.96  2006/08/25 16:56:30  kharlov
  80  * Compliance with Effective C++
  81  *
  82  * Revision 1.95  2006/08/11 12:36:26  cvetan
  83  * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
  84  *
  85  * Revision 1.94  2006/08/07 12:27:49  hristov
  86  * Removing obsolete code which affected the event numbering scheme
  87  *
  88  * Revision 1.93  2006/08/01 12:20:17  cvetan
  89  * 1. Adding a possibility to read and reconstruct an old rcu formatted raw data. This is controlled by an option of AliReconstruction and AliPHOSReconstructor. 2. In case of raw data processing (without galice.root) create the default AliPHOSGeometry object. Most likely this should be moved to the CDB
  90  *
  91  * Revision 1.92  2006/04/29 20:26:46  hristov
  92  * Separate EMC and CPV calibration (Yu.Kharlov)
  93  *
  94  * Revision 1.91  2006/04/22 10:30:17  hristov
  95  * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
  96  *
  97  * Revision 1.90  2006/04/11 15:22:59  hristov
  98  * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
  99  *
 100  * Revision 1.89  2006/03/13 14:05:42  kharlov
 101  * Calibration objects for EMC and CPV
 102  *
 103  * Revision 1.88  2006/01/11 08:54:52  hristov
 104  * Additional protection in case no calibration entry was found
 105  *
 106  * Revision 1.87  2005/11/22 08:46:43  kharlov
 107  * Updated to new CDB (Boris Polichtchouk)
 108  *
 109  * Revision 1.86  2005/11/14 21:52:43  hristov
 110  * Coding conventions
 111  *
 112  * Revision 1.85  2005/09/27 16:08:08  hristov
 113  * New version of CDB storage framework (A.Colla)
 114  *
 115  * Revision 1.84  2005/09/21 10:02:47  kharlov
 116  * Reading calibration from CDB (Boris Polichtchouk)
 117  *
 118  * Revision 1.82  2005/09/02 15:43:13  kharlov
 119  * Add comments in GetCalibrationParameters and Calibrate
 120  *
 121  * Revision 1.81  2005/09/02 14:32:07  kharlov
 122  * Calibration of raw data
 123  *
 124  * Revision 1.80  2005/08/24 15:31:36  kharlov
 125  * Setting raw digits flag
 126  *
 127  * Revision 1.79  2005/07/25 15:53:53  kharlov
 128  * Read raw data
 129  *
 130  * Revision 1.78  2005/05/28 14:19:04  schutz
 131  * Compilation warnings fixed by T.P.
 132  *
 133  */
 134
 135 //*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
 136 //////////////////////////////////////////////////////////////////////////////
 137 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and
 138 //  unfolds the clusters having several local maxima.
 139 //  Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
 140 //  PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
 141 //  parameters including input digits branch title, thresholds etc.)
 142 //  This TTask is normally called from Reconstructioner, but can as well be used in
 143 //  standalone mode.
 144 // Use Case:
 145 //  root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1(<pointer_to_phos_geometry_onject>)
 146 //  root [1] cl->Digits2Clusters(digitsTree,clusterTree)
 147 //               //finds RecPoints in the current event
 148 //  root [2] cl->SetDigitsBranch("digits2")
 149 //               //sets another title for Digitis (input) branch
 150 //  root [3] cl->SetRecPointsBranch("recp2")
 151 //               //sets another title four output branches
 152 //  root [4] cl->SetEmcLocalMaxCut(0.03)
 153 //               //set clusterization parameters
 154
 155 // --- ROOT system ---
 156
 157 #include "TMath.h"
 158 #include "TMinuit.h"
 159 #include "TTree.h"
 160 #include "TBenchmark.h"
 161 #include "TClonesArray.h"
 162
 163 // --- Standard library ---
 164
 165 // --- AliRoot header files ---
 166 #include "AliLog.h"
 167 #include "AliConfig.h"
 168 #include "AliPHOSGeometry.h"
 169 #include "AliPHOSClusterizerv1.h"
 170 #include "AliPHOSEmcRecPoint.h"
 171 #include "AliPHOSCpvRecPoint.h"
 172 #include "AliPHOSDigit.h"
 173 #include "AliPHOSDigitizer.h"
 174 #include "AliPHOSCalibrationDB.h"
 175 #include "AliCDBManager.h"
 176 #include "AliCDBStorage.h"
 177 #include "AliCDBEntry.h"
 178 #include "AliPHOSRecoParam.h"
 179 #include "AliPHOSCalibData.h"
 180 #include "AliPHOSReconstructor.h"
 181
 182 ClassImp(AliPHOSClusterizerv1)
 183
 184 //____________________________________________________________________________
 185 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
 186   AliPHOSClusterizer(),
 187   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 188   fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
 189   fADCchanelEmc(0),         fADCpedestalEmc(0),
 190   fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
 191   fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
 192   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 193   fW0CPV(0),                  fEmcTimeGate(0),
 194   fIsOldRCUFormat(0)
 195 {
 196   // default ctor (to be used mainly by Streamer)
 197
 198   InitParameters() ;
 199   fDefaultInit = kTRUE ;
 200 }
 201
 202 //____________________________________________________________________________
 203 AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
 204   AliPHOSClusterizer(geom),
 205   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 206   fWrite(0),                  fNumberOfEmcClusters(0),  fNumberOfCpvClusters(0),
 207   fADCchanelEmc(0),         fADCpedestalEmc(0),
 208   fADCchanelCpv(0),           fADCpedestalCpv(0),       fEmcClusteringThreshold(0),
 209   fCpvClusteringThreshold(0), fEmcMinE(0),              fCpvMinE(0),
 210   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 211   fW0CPV(0),                  fEmcTimeGate(0),
 212   fIsOldRCUFormat(0)
 213 {
 214   // ctor with the indication of the file where header Tree and digits Tree are stored
 215
 216   InitParameters() ;
 217   Init() ;
 218   fDefaultInit = kFALSE ;
 219 }
 220
 221 //____________________________________________________________________________
 222   AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 223 {
 224   // dtor
 225
 226 }
 227
 228 //____________________________________________________________________________
 229 Float_t  AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
 230 {
 231   // Convert EMC measured amplitude into real energy.
 232   // Calibration parameters are taken from calibration data base for raw data,
 233   // or from digitizer parameters for simulated data.
 234
 235   if(fgCalibData){
 236     Int_t relId[4];
 237     fGeom->AbsToRelNumbering(absId,relId) ;
 238     Int_t   module = relId[0];
 239     Int_t   column = relId[3];
 240     Int_t   row    = relId[2];
 241     if(absId <= fEmcCrystals) { // this is EMC
 242       fADCchanelEmc   = fgCalibData->GetADCchannelEmc (module,column,row);
 243       return amp*fADCchanelEmc ;
 244     }
 245   }
 246   else{ //simulation
 247     if(absId <= fEmcCrystals) // this is EMC
 248       return fADCpedestalEmc + amp*fADCchanelEmc ;
 249   }
 250   return 0;
 251 }
 252
 253 //____________________________________________________________________________
 254 Float_t  AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
 255 {
 256   // Convert digitized CPV amplitude into charge.
 257   // Calibration parameters are taken from calibration data base for raw data,
 258   // or from digitizer parameters for simulated data.
 259
 260   if(fgCalibData){
 261     Int_t relId[4];
 262     fGeom->AbsToRelNumbering(absId,relId) ;
 263     Int_t   module = relId[0];
 264     Int_t   column = relId[3];
 265     Int_t   row    = relId[2];
 266     if(absId > fEmcCrystals) { // this is CPV
 267       fADCchanelCpv   = fgCalibData->GetADCchannelCpv (module,column,row);
 268       fADCpedestalCpv = fgCalibData->GetADCpedestalCpv(module,column,row);
 269       return fADCpedestalCpv + amp*fADCchanelCpv ;
 270     }
 271   }
 272   else{ //simulation
 273     if(absId > fEmcCrystals) // this is CPV
 274       return fADCpedestalCpv+ amp*fADCchanelCpv ;
 275   }
 276   return 0;
 277 }
 278
 279 //____________________________________________________________________________
 280 void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
 281 {
 282   // Steering method to perform clusterization for one event
 283   // The input is the tree with digits.
 284   // The output is the tree with clusters.
 285
 286   if(strstr(option,"tim"))
 287     gBenchmark->Start("PHOSClusterizer");
 288
 289   if(strstr(option,"print")) {
 290     Print() ;
 291     return ;
 292   }
 293
 294   MakeClusters() ;
 295
 296   AliDebug(2,Form(" ---- Printing clusters (%d)\n",
 297                   fEMCRecPoints->GetEntries()));
 298   if(AliLog::GetGlobalDebugLevel()>1)
 299     fEMCRecPoints->Print();
 300
 301   if(fToUnfold)
 302     MakeUnfolding();
 303
 304   WriteRecPoints();
 305
 306   if(strstr(option,"deb"))
 307     PrintRecPoints(option) ;
 308
 309   if(strstr(option,"tim")){
 310     gBenchmark->Stop("PHOSClusterizer");
 311     AliInfo(Form("took %f seconds for Clusterizing\n",
 312                  gBenchmark->GetCpuTime("PHOSClusterizer")));
 313   }
 314 }
 315
 316 //____________________________________________________________________________
 317 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
 318                                     Int_t nPar, Float_t * fitparameters) const
 319 {
 320   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 321   // The initial values for fitting procedure are set equal to the positions of local maxima.
 322   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 323
 324
 325   gMinuit->mncler();                     // Reset Minuit's list of paramters
 326   gMinuit->SetPrintLevel(-1) ;           // No Printout
 327   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
 328                                          // To set the address of the minimization function
 329
 330   TList * toMinuit = new TList();
 331   toMinuit->AddAt(emcRP,0) ;
 332   toMinuit->AddAt(fDigitsArr,1) ;
 333   toMinuit->AddAt(fGeom,2) ;
 334
 335   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 336
 337   // filling initial values for fit parameters
 338   AliPHOSDigit * digit ;
 339
 340   Int_t ierflg  = 0;
 341   Int_t index   = 0 ;
 342   Int_t nDigits = (Int_t) nPar / 3 ;
 343
 344   Int_t iDigit ;
 345
 346   for(iDigit = 0; iDigit < nDigits; iDigit++){
 347     digit = maxAt[iDigit];
 348
 349     Int_t relid[4] ;
 350     Float_t x = 0.;
 351     Float_t z = 0.;
 352     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 353     fGeom->RelPosInModule(relid, x, z) ;
 354
 355     Float_t energy = maxAtEnergy[iDigit] ;
 356
 357     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 358     index++ ;
 359     if(ierflg != 0){
 360       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
 361       return kFALSE;
 362     }
 363     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 364     index++ ;
 365     if(ierflg != 0){
 366        Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
 367       return kFALSE;
 368     }
 369     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 370     index++ ;
 371     if(ierflg != 0){
 372       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
 373       return kFALSE;
 374     }
 375   }
 376
 377   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 378                       //  depends on it.
 379   Double_t p1 = 1.0 ;
 380   Double_t p2 = 0.0 ;
 381
 382   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 383   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 384   gMinuit->SetMaxIterations(5);
 385   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 386
 387   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 388
 389   if(ierflg == 4){  // Minimum not found
 390     Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
 391     return kFALSE ;
 392   }
 393   for(index = 0; index < nPar; index++){
 394     Double_t err ;
 395     Double_t val ;
 396     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 397     fitparameters[index] = val ;
 398    }
 399
 400   delete toMinuit ;
 401   return kTRUE;
 402
 403 }
 404
 405 //____________________________________________________________________________
 406 void AliPHOSClusterizerv1::GetCalibrationParameters()
 407 {
 408   // Set calibration parameters:
 409   // if calibration database exists, they are read from database,
 410   // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
 411   //
 412   // It is a user responsilibity to open CDB before reconstruction, for example:
 413   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 414
 415   if (!fgCalibData)
 416     fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
 417   if (fgCalibData->GetCalibDataEmc() == 0)
 418     AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
 419   if (fgCalibData->GetCalibDataCpv() == 0)
 420     AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
 421
 422 }
 423
 424 //____________________________________________________________________________
 425 void AliPHOSClusterizerv1::Init()
 426 {
 427   // Make all memory allocations which can not be done in default constructor.
 428   // Attach the Clusterizer task to the list of PHOS tasks
 429
 430   fEmcCrystals = fGeom->GetNModules() *  fGeom->GetNCristalsInModule() ;
 431
 432   if(!gMinuit)
 433     gMinuit = new TMinuit(100);
 434
 435   GetCalibrationParameters() ;
 436
 437 }
 438
 439 //____________________________________________________________________________
 440 void AliPHOSClusterizerv1::InitParameters()
 441 {
 442
 443   fNumberOfCpvClusters     = 0 ;
 444   fNumberOfEmcClusters     = 0 ;
 445
 446   const AliPHOSRecoParam* parEmc = AliPHOSReconstructor::GetRecoParamEmc();
 447   if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!");
 448
 449   const AliPHOSRecoParam* parCpv = AliPHOSReconstructor::GetRecoParamCpv();
 450   if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!");
 451
 452   fCpvClusteringThreshold  = parCpv->GetClusteringThreshold();
 453   fEmcClusteringThreshold  = parEmc->GetClusteringThreshold();
 454
 455   fEmcLocMaxCut            = parEmc->GetLocalMaxCut();
 456   fCpvLocMaxCut            = parCpv->GetLocalMaxCut();
 457
 458   fEmcMinE                 = parEmc->GetMinE();
 459   fCpvMinE                 = parCpv->GetMinE();
 460
 461   fW0                      = parEmc->GetLogWeight();
 462   fW0CPV                   = parCpv->GetLogWeight();
 463
 464   fEmcTimeGate             = 1.e-8 ;
 465
 466   fToUnfold                = kTRUE ;
 467
 468   fWrite                   = kTRUE ;
 469
 470   fIsOldRCUFormat          = kFALSE;
 471 }
 472
 473 //____________________________________________________________________________
 474 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
 475 {
 476   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 477   //                                       = 1 are neighbour
 478   //                                       = 2 are not neighbour but do not continue searching
 479   // neighbours are defined as digits having at least a common vertex
 480   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 481   //                                      which is compared to a digit (d2)  not yet in a cluster
 482
 483   Int_t rv = 0 ;
 484
 485   Int_t relid1[4] ;
 486   fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
 487
 488   Int_t relid2[4] ;
 489   fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
 490
 491   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
 492     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 493     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 494
 495     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
 496       if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
 497       rv = 1 ;
 498     }
 499     else {
 500       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 501         rv = 2; //  Difference in row numbers is too large to look further
 502     }
 503
 504   }
 505   else {
 506
 507     if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
 508       rv=2 ;
 509
 510   }
 511
 512   return rv ;
 513 }
 514 //____________________________________________________________________________
 515 void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
 516 {
 517   // Remove digits with amplitudes below threshold
 518
 519   for(Int_t i=0; i<digits->GetEntriesFast(); i++){
 520     AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
 521     if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
 522          (IsInCpv(digit) && CalibrateCPV(digit->GetAmp()   ,digit->GetId()) < fCpvMinE) )
 523       digits->RemoveAt(i) ;
 524   }
 525   digits->Compress() ;
 526   for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
 527     AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
 528     digit->SetIndexInList(i) ;
 529   }
 530 }
 531
 532 //____________________________________________________________________________
 533 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 534 {
 535   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
 536
 537   Bool_t rv = kFALSE ;
 538
 539   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 540
 541   if(digit->GetId() <= nEMC )   rv = kTRUE;
 542
 543   return rv ;
 544 }
 545
 546 //____________________________________________________________________________
 547 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 548 {
 549   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
 550
 551   Bool_t rv = kFALSE ;
 552
 553   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 554
 555   if(digit->GetId() > nEMC )   rv = kTRUE;
 556
 557   return rv ;
 558 }
 559
 560 //____________________________________________________________________________
 561 void AliPHOSClusterizerv1::WriteRecPoints()
 562 {
 563
 564   // Creates new branches with given title
 565   // fills and writes into TreeR.
 566
 567   Int_t index ;
 568   //Evaluate position, dispersion and other RecPoint properties..
 569   Int_t nEmc = fEMCRecPoints->GetEntriesFast();
 570   for(index = 0; index < nEmc; index++){
 571     AliPHOSEmcRecPoint * rp =
 572       dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
 573     rp->Purify(fEmcMinE) ;
 574     if(rp->GetMultiplicity()==0){
 575       fEMCRecPoints->RemoveAt(index) ;
 576       delete rp ;
 577       continue;
 578     }
 579
 580     // No vertex is available now, calculate corrections in PID
 581     rp->EvalAll(fW0,fDigitsArr) ;
 582     TVector3 fakeVtx(0.,0.,0.) ;
 583     rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
 584     rp->EvalLocal2TrackingCSTransform();
 585   }
 586   fEMCRecPoints->Compress() ;
 587   fEMCRecPoints->Sort() ;
 588   //  fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
 589   for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
 590     dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
 591   }
 592
 593   //For each rec.point set the distance to the nearest bad crystal (BVP)
 594   SetDistancesToBadChannels();
 595
 596   //Now the same for CPV
 597   for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
 598     AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
 599     rp->EvalAll(fW0CPV,fDigitsArr) ;
 600     rp->EvalLocal2TrackingCSTransform();
 601   }
 602   fCPVRecPoints->Sort() ;
 603
 604   for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
 605     dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
 606
 607   fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
 608
 609   if(fWrite){ //We write TreeR
 610     fTreeR->Fill();
 611   }
 612 }
 613
 614 //____________________________________________________________________________
 615 void AliPHOSClusterizerv1::MakeClusters()
 616 {
 617   // Steering method to construct the clusters stored in a list of Reconstructed Points
 618   // A cluster is defined as a list of neighbour digits
 619
 620   //Remove digits below threshold
 621   CleanDigits(fDigitsArr) ;
 622
 623   TClonesArray * digitsC =  static_cast<TClonesArray*>( fDigitsArr->Clone() ) ;
 624
 625   // Clusterization starts
 626
 627   TIter nextdigit(digitsC) ;
 628   AliPHOSDigit * digit ;
 629   Bool_t notremoved = kTRUE ;
 630
 631   while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
 632
 633
 634     AliPHOSRecPoint * clu = 0 ;
 635
 636     TArrayI clusterdigitslist(1500) ;
 637     Int_t index ;
 638
 639     if (( IsInEmc (digit) &&
 640           CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
 641         ( IsInCpv (digit) &&
 642           CalibrateCPV(digit->GetAmp()   ,digit->GetId()) > fCpvClusteringThreshold ) ) {
 643       Int_t iDigitInCluster = 0 ;
 644
 645       if  ( IsInEmc(digit) ) {
 646         // start a new EMC RecPoint
 647         if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 648           fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
 649
 650         fEMCRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
 651         clu = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 652         fNumberOfEmcClusters++ ;
 653         clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
 654         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
 655         iDigitInCluster++ ;
 656         digitsC->Remove(digit) ;
 657
 658       } else {
 659
 660         // start a new CPV cluster
 661         if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 662           fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
 663
 664         fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
 665
 666         clu =  dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 667         fNumberOfCpvClusters++ ;
 668         clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
 669         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;
 670         iDigitInCluster++ ;
 671         digitsC->Remove(digit) ;
 672         nextdigit.Reset() ;
 673
 674         // Here we remove remaining EMC digits, which cannot make a cluster
 675
 676         if( notremoved ) {
 677           while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
 678             if( IsInEmc(digit) )
 679               digitsC->Remove(digit) ;
 680             else
 681               break ;
 682           }
 683           notremoved = kFALSE ;
 684         }
 685
 686       } // else
 687
 688       nextdigit.Reset() ;
 689
 690       AliPHOSDigit * digitN ;
 691       index = 0 ;
 692       while (index < iDigitInCluster){ // scan over digits already in cluster
 693         digit =  dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) )  ;
 694         index++ ;
 695         while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
 696           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 697           switch (ineb ) {
 698           case 0 :   // not a neighbour
 699             break ;
 700           case 1 :   // are neighbours
 701             if (IsInEmc (digitN))
 702               clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
 703             else
 704               clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp()   , digitN->GetId() ) );
 705
 706             clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
 707             iDigitInCluster++ ;
 708             digitsC->Remove(digitN) ;
 709             break ;
 710           case 2 :   // too far from each other
 711             goto endofloop;
 712           } // switch
 713
 714         } // while digitN
 715
 716       endofloop: ;
 717         nextdigit.Reset() ;
 718
 719       } // loop over cluster
 720
 721     } // energy theshold
 722
 723
 724   } // while digit
 725
 726   delete digitsC ;
 727
 728 }
 729
 730 //____________________________________________________________________________
 731 void AliPHOSClusterizerv1::MakeUnfolding()
 732 {
 733   // Unfolds clusters using the shape of an ElectroMagnetic shower
 734   // Performs unfolding of all EMC/CPV clusters
 735
 736   // Unfold first EMC clusters
 737   if(fNumberOfEmcClusters > 0){
 738
 739     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 740
 741     Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
 742     Int_t index ;
 743     for(index = 0 ; index < numberofNotUnfolded ; index++){
 744
 745       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
 746       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
 747         break ;
 748
 749       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 750       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 751       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 752       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
 753
 754       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 755         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 756         fEMCRecPoints->Remove(emcRecPoint);
 757         fEMCRecPoints->Compress() ;
 758         index-- ;
 759         fNumberOfEmcClusters -- ;
 760         numberofNotUnfolded-- ;
 761       }
 762       else{
 763         emcRecPoint->SetNExMax(1) ; //Only one local maximum
 764       }
 765
 766       delete[] maxAt ;
 767       delete[] maxAtEnergy ;
 768     }
 769   }
 770   // Unfolding of EMC clusters finished
 771
 772
 773   // Unfold now CPV clusters
 774   if(fNumberOfCpvClusters > 0){
 775
 776     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 777
 778     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
 779     Int_t index ;
 780     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
 781
 782       AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
 783
 784       if(recPoint->GetPHOSMod()> nModulesToUnfold)
 785         break ;
 786
 787       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
 788
 789       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 790       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 791       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 792       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
 793
 794       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 795         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 796         fCPVRecPoints->Remove(emcRecPoint);
 797         fCPVRecPoints->Compress() ;
 798         index-- ;
 799         numberofCpvNotUnfolded-- ;
 800         fNumberOfCpvClusters-- ;
 801       }
 802
 803       delete[] maxAt ;
 804       delete[] maxAtEnergy ;
 805     }
 806   }
 807   //Unfolding of Cpv clusters finished
 808
 809 }
 810
 811 //____________________________________________________________________________
 812 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
 813 {
 814   // Shape of the shower (see PHOS TDR)
 815   // If you change this function, change also the gradient evaluation in ChiSquare()
 816
 817   //for the moment we neglect dependence on the incident angle.
 818
 819   Double_t r2    = x*x + z*z ;
 820   Double_t r4    = r2*r2 ;
 821   Double_t r295  = TMath::Power(r2, 2.95/2.) ;
 822   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 823   return shape ;
 824 }
 825
 826 //____________________________________________________________________________
 827 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 828                                           Int_t nMax,
 829                                           AliPHOSDigit ** maxAt,
 830                                           Float_t * maxAtEnergy)
 831 {
 832   // Performs the unfolding of a cluster with nMax overlapping showers
 833
 834   Int_t nPar = 3 * nMax ;
 835   Float_t * fitparameters = new Float_t[nPar] ;
 836
 837   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 838   if( !rv ) {
 839     // Fit failed, return and remove cluster
 840     iniEmc->SetNExMax(-1) ;
 841     delete[] fitparameters ;
 842     return ;
 843   }
 844
 845   // create ufolded rec points and fill them with new energy lists
 846   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
 847   // and later correct this number in acordance with actual energy deposition
 848
 849   Int_t nDigits = iniEmc->GetMultiplicity() ;
 850   Float_t * efit = new Float_t[nDigits] ;
 851   Float_t xDigit=0.,zDigit=0. ;
 852   Float_t xpar=0.,zpar=0.,epar=0.  ;
 853   Int_t relid[4] ;
 854   AliPHOSDigit * digit = 0 ;
 855   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 856
 857   TVector3 vIncid ;
 858
 859   Int_t iparam ;
 860   Int_t iDigit ;
 861   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 862     digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
 863     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 864     fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 865     efit[iDigit] = 0;
 866
 867     iparam = 0 ;
 868     while(iparam < nPar ){
 869       xpar = fitparameters[iparam] ;
 870       zpar = fitparameters[iparam+1] ;
 871       epar = fitparameters[iparam+2] ;
 872       iparam += 3 ;
 873 //      fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
 874 //      efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
 875       efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 876     }
 877   }
 878
 879
 880   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
 881   // so that energy deposited in each cell is distributed betwin new clusters proportionally
 882   // to its contribution to efit
 883
 884   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 885   Float_t ratio ;
 886
 887   iparam = 0 ;
 888   while(iparam < nPar ){
 889     xpar = fitparameters[iparam] ;
 890     zpar = fitparameters[iparam+1] ;
 891     epar = fitparameters[iparam+2] ;
 892     iparam += 3 ;
 893 //    fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
 894
 895     AliPHOSEmcRecPoint * emcRP = 0 ;
 896
 897     if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
 898
 899       if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 900         fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
 901
 902       (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
 903       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 904       fNumberOfEmcClusters++ ;
 905       emcRP->SetNExMax((Int_t)nPar/3) ;
 906     }
 907     else{//create new entries in fCPVRecPoints
 908       if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 909         fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
 910
 911       (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
 912       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 913       fNumberOfCpvClusters++ ;
 914     }
 915
 916     Float_t eDigit ;
 917     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 918       digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
 919       fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 920       fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 921 //      ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
 922       ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
 923       eDigit = emcEnergies[iDigit] * ratio ;
 924       emcRP->AddDigit( *digit, eDigit ) ;
 925     }
 926   }
 927
 928   delete[] fitparameters ;
 929   delete[] efit ;
 930
 931 }
 932
 933 //_____________________________________________________________________________
 934 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 935 {
 936   // Calculates the Chi square for the cluster unfolding minimization
 937   // Number of parameters, Gradient, Chi squared, parameters, what to do
 938
 939   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
 940
 941   AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
 942   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
 943   // A bit buggy way to get an access to the geometry
 944   // To be revised!
 945   AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
 946
 947 //  TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ;  //Vertex position
 948
 949   //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
 950
 951   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 952
 953   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
 954
 955   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
 956
 957 //  TVector3 vInc ;
 958   fret = 0. ;
 959   Int_t iparam ;
 960
 961   if(iflag == 2)
 962     for(iparam = 0 ; iparam < nPar ; iparam++)
 963       Grad[iparam] = 0 ; // Will evaluate gradient
 964
 965   Double_t efit ;
 966
 967   AliPHOSDigit * digit ;
 968   Int_t iDigit ;
 969
 970   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 971
 972     digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
 973
 974     Int_t relid[4] ;
 975     Float_t xDigit ;
 976     Float_t zDigit ;
 977
 978     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 979
 980     geom->RelPosInModule(relid, xDigit, zDigit) ;
 981
 982      if(iflag == 2){  // calculate gradient
 983        Int_t iParam = 0 ;
 984        efit = 0 ;
 985        while(iParam < nPar ){
 986          Double_t dx = (xDigit - x[iParam]) ;
 987          iParam++ ;
 988          Double_t dz = (zDigit - x[iParam]) ;
 989          iParam++ ;
 990 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
 991 //         efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
 992          efit += x[iParam] * ShowerShape(dx,dz) ;
 993          iParam++ ;
 994        }
 995        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
 996        iParam = 0 ;
 997        while(iParam < nPar ){
 998          Double_t xpar = x[iParam] ;
 999          Double_t zpar = x[iParam+1] ;
1000          Double_t epar = x[iParam+2] ;
1001 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
1002          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1003 //         Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1004          Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1005 //DP: No incident angle dependence in gradient yet!!!!!!
1006          Double_t r4 = dr*dr*dr*dr ;
1007          Double_t r295 = TMath::Power(dr,2.95) ;
1008          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1009                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1010
1011          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
1012          iParam++ ;
1013          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
1014          iParam++ ;
1015          Grad[iParam] += shape ;                                  // Derivative over energy
1016          iParam++ ;
1017        }
1018      }
1019      efit = 0;
1020      iparam = 0 ;
1021
1022      while(iparam < nPar ){
1023        Double_t xpar = x[iparam] ;
1024        Double_t zpar = x[iparam+1] ;
1025        Double_t epar = x[iparam+2] ;
1026        iparam += 3 ;
1027 //       fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
1028 //       efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1029        efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1030      }
1031
1032      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
1033      // Here we assume, that sigma = sqrt(E)
1034   }
1035
1036 }
1037
1038 //____________________________________________________________________________
1039 void AliPHOSClusterizerv1::Print(const Option_t *)const
1040 {
1041   // Print clusterizer parameters
1042
1043   TString message ;
1044   TString taskName(GetName()) ;
1045   taskName.ReplaceAll(Version(), "") ;
1046
1047   if( strcmp(GetName(), "") !=0 ) {
1048     // Print parameters
1049     message  = "\n--------------- %s %s -----------\n" ;
1050     message += "Clusterizing digits from the file: %s\n" ;
1051     message += "                           Branch: %s\n" ;
1052     message += "                       EMC Clustering threshold = %f\n" ;
1053     message += "                       EMC Local Maximum cut    = %f\n" ;
1054     message += "                       EMC Logarothmic weight   = %f\n" ;
1055     message += "                       CPV Clustering threshold = %f\n" ;
1056     message += "                       CPV Local Maximum cut    = %f\n" ;
1057     message += "                       CPV Logarothmic weight   = %f\n" ;
1058     if(fToUnfold)
1059       message += " Unfolding on\n" ;
1060     else
1061       message += " Unfolding off\n" ;
1062
1063     message += "------------------------------------------------------------------" ;
1064   }
1065   else
1066     message = " AliPHOSClusterizerv1 not initialized " ;
1067
1068   AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
1069        taskName.Data(),
1070        GetTitle(),
1071        taskName.Data(),
1072        GetName(),
1073        fEmcClusteringThreshold,
1074        fEmcLocMaxCut,
1075        fW0,
1076        fCpvClusteringThreshold,
1077        fCpvLocMaxCut,
1078        fW0CPV )) ;
1079 }
1080 //____________________________________________________________________________
1081 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1082 {
1083   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1084
1085   AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
1086                fEMCRecPoints->GetEntriesFast(),
1087                fCPVRecPoints->GetEntriesFast() ))  ;
1088
1089   if(strstr(option,"all")) {
1090     printf("\n EMC clusters \n") ;
1091     printf("Index    Ene(MeV) Multi Module     X    Y   Z    Lambdas_1  Lambda_2  # of prim  Primaries list\n") ;
1092     Int_t index ;
1093     for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1094       AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1095       TVector3  locpos;
1096       rp->GetLocalPosition(locpos);
1097       Float_t lambda[2];
1098       rp->GetElipsAxis(lambda);
1099       Int_t * primaries;
1100       Int_t nprimaries;
1101       primaries = rp->GetPrimaries(nprimaries);
1102       printf("\n%6d  %8.2f  %3d     %2d     %4.1f   %4.1f   %4.1f   %4f  %4f    %2d     : ",
1103               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1104               locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1105
1106       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1107         printf("%d ", primaries[iprimary] ) ;
1108       }
1109       printf("\n") ;
1110     }
1111
1112     //Now plot CPV recPoints
1113     printf("\n CPV clusters \n") ;
1114     printf("Index    Ene(MeV) Module     X     Y    Z  \n") ;
1115     for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1116       AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1117
1118       TVector3  locpos;
1119       rp->GetLocalPosition(locpos);
1120
1121       printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f \n",
1122              rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1123              locpos.X(), locpos.Y(), locpos.Z()) ;
1124     }
1125   }
1126 }
1127
1128
1129 //____________________________________________________________________________
1130 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1131 {
1132   //For each EMC rec. point set the distance to the nearest bad crystal.
1133   //Author: Boris Polichtchouk
1134
1135   if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1136   AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
1137
1138   Int_t badIds[8000];
1139   fgCalibData->EmcBadChannelIds(badIds);
1140
1141   AliPHOSEmcRecPoint* rp;
1142
1143   TMatrixF gmat;
1144   TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1145   TVector3 gposBadChannel; // global position of bad crystal
1146   TVector3 dR;
1147
1148   Float_t dist,minDist;
1149
1150   for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1151     rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1152     minDist = 1.e+07;
1153
1154     for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1155       rp->GetGlobalPosition(gposRecPoint,gmat);
1156       fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1157       AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1158                       gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1159                       gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1160       dR = gposBadChannel-gposRecPoint;
1161       dist = dR.Mag();
1162       if(dist<minDist) minDist = dist;
1163     }
1164
1165     rp->SetDistanceToBadCrystal(minDist);
1166   }
1167
1168 }