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