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