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