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