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