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: 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
 209 //____________________________________________________________________________
 210 AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
 211   AliPHOSClusterizer(geom),
 212   fDefaultInit(0),            fEmcCrystals(0),          fToUnfold(0),
 213   fWrite(0),
 214   fNumberOfEmcClusters(0),    fNumberOfCpvClusters(0),
 215   fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
 216   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
 217   fW0CPV(0),
 218   fTimeGateLowAmp(0.),        fTimeGateLow(0.),         fTimeGateHigh(0.),
 219   fEcoreRadius(0)
 220 {
 221   // ctor with the indication of the file where header Tree and digits Tree are stored
 222
 223   Init() ;
 224   fDefaultInit = kFALSE ;
 225 }
 226
 227 //____________________________________________________________________________
 228   AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
 229 {
 230   // dtor
 231
 232 }
 233 //____________________________________________________________________________
 234 void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
 235 {
 236   // Steering method to perform clusterization for one event
 237   // The input is the tree with digits.
 238   // The output is the tree with clusters.
 239
 240   if(strstr(option,"tim"))
 241     gBenchmark->Start("PHOSClusterizer");
 242
 243   if(strstr(option,"print")) {
 244     Print() ;
 245     return ;
 246   }
 247
 248   MakeClusters() ;
 249
 250   AliDebug(2,Form(" ---- Printing clusters (%d)\n",
 251                   fEMCRecPoints->GetEntries()));
 252   if(AliLog::GetGlobalDebugLevel()>1)
 253     fEMCRecPoints->Print();
 254
 255   if(fToUnfold)
 256     MakeUnfolding();
 257
 258   WriteRecPoints();
 259
 260   if(strstr(option,"deb"))
 261     PrintRecPoints(option) ;
 262
 263   if(strstr(option,"tim")){
 264     gBenchmark->Stop("PHOSClusterizer");
 265     AliInfo(Form("took %f seconds for Clusterizing\n",
 266                  gBenchmark->GetCpuTime("PHOSClusterizer")));
 267   }
 268   fEMCRecPoints->Clear("C");
 269   fCPVRecPoints->Clear("C");
 270 }
 271
 272 //____________________________________________________________________________
 273 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
 274                                     Int_t nPar, Float_t * fitparameters) const
 275 {
 276   // Calls TMinuit to fit the energy distribution of a cluster with several maxima
 277   // The initial values for fitting procedure are set equal to the positions of local maxima.
 278   // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
 279
 280
 281   if(!gMinuit) //it was deleted by someone else
 282     gMinuit = new TMinuit(100) ;
 283   gMinuit->mncler();                     // Reset Minuit's list of paramters
 284   gMinuit->SetPrintLevel(-1) ;           // No Printout
 285   gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
 286                                          // To set the address of the minimization function
 287
 288   TList * toMinuit = new TList();
 289   toMinuit->AddAt(emcRP,0) ;
 290   toMinuit->AddAt(fDigitsArr,1) ;
 291   toMinuit->AddAt(fGeom,2) ;
 292
 293   gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
 294
 295   // filling initial values for fit parameters
 296   AliPHOSDigit * digit ;
 297
 298   Int_t ierflg  = 0;
 299   Int_t index   = 0 ;
 300   Int_t nDigits = (Int_t) nPar / 3 ;
 301
 302   Int_t iDigit ;
 303
 304   for(iDigit = 0; iDigit < nDigits; iDigit++){
 305     digit = maxAt[iDigit];
 306
 307     Int_t relid[4] ;
 308     Float_t x = 0.;
 309     Float_t z = 0.;
 310     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 311     fGeom->RelPosInModule(relid, x, z) ;
 312
 313     Float_t energy = maxAtEnergy[iDigit] ;
 314
 315     gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
 316     index++ ;
 317     if(ierflg != 0){
 318       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
 319       return kFALSE;
 320     }
 321     gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
 322     index++ ;
 323     if(ierflg != 0){
 324        Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
 325       return kFALSE;
 326     }
 327     gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
 328     index++ ;
 329     if(ierflg != 0){
 330       Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
 331       return kFALSE;
 332     }
 333   }
 334
 335   Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
 336                       //  depends on it.
 337   Double_t p1 = 1.0 ;
 338   Double_t p2 = 0.0 ;
 339
 340   gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
 341   gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
 342   gMinuit->SetMaxIterations(5);
 343   gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
 344
 345   gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
 346
 347   if(ierflg == 4){  // Minimum not found
 348     Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
 349     return kFALSE ;
 350   }
 351   for(index = 0; index < nPar; index++){
 352     Double_t err ;
 353     Double_t val ;
 354     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
 355     fitparameters[index] = val ;
 356    }
 357
 358   delete toMinuit ;
 359   return kTRUE;
 360
 361 }
 362
 363
 364 //____________________________________________________________________________
 365 void AliPHOSClusterizerv1::Init()
 366 {
 367   // Make all memory allocations which can not be done in default constructor.
 368   // Attach the Clusterizer task to the list of PHOS tasks
 369
 370   fEmcCrystals = fGeom->GetNModules() *  fGeom->GetNCristalsInModule() ;
 371
 372   if(!gMinuit)
 373     gMinuit = new TMinuit(100);
 374
 375   if (!fgCalibData)
 376     fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
 377   if (fgCalibData->GetCalibDataEmc() == 0)
 378     AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
 379   if (fgCalibData->GetCalibDataCpv() == 0)
 380     AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
 381
 382 }
 383
 384 //____________________________________________________________________________
 385 void AliPHOSClusterizerv1::InitParameters()
 386 {
 387
 388   fNumberOfCpvClusters     = 0 ;
 389   fNumberOfEmcClusters     = 0 ;
 390
 391   const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
 392   if(!recoParam) AliFatal("Reconstruction parameters are not set!");
 393
 394   recoParam->Print();
 395
 396   fEmcClusteringThreshold  = recoParam->GetEMCClusteringThreshold();
 397   fCpvClusteringThreshold  = recoParam->GetCPVClusteringThreshold();
 398
 399   fEmcLocMaxCut            = recoParam->GetEMCLocalMaxCut();
 400   fCpvLocMaxCut            = recoParam->GetCPVLocalMaxCut();
 401
 402   fW0                      = recoParam->GetEMCLogWeight();
 403   fW0CPV                   = recoParam->GetCPVLogWeight();
 404
 405   fTimeGateLowAmp          = recoParam->GetTimeGateAmpThresh() ;
 406   fTimeGateLow             = recoParam->GetTimeGateLow() ;
 407   fTimeGateHigh            = recoParam->GetTimeGateHigh() ;
 408
 409   fEcoreRadius             = recoParam->GetEMCEcoreRadius();
 410
 411   fToUnfold                = recoParam->EMCToUnfold() ;
 412
 413   fWrite                   = kTRUE ;
 414 }
 415
 416 //____________________________________________________________________________
 417 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
 418 {
 419   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
 420   //                                       = 1 are neighbour
 421   //                                       = 2 are not neighbour but do not continue searching
 422   //                                       =-1 are not neighbour, continue searching, but do not look before d2 next time
 423   // neighbours are defined as digits having at least a common vertex
 424   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 425   //                                      which is compared to a digit (d2)  not yet in a cluster
 426
 427   Int_t relid1[4] ;
 428   fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
 429
 430   Int_t relid2[4] ;
 431   fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
 432
 433   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
 434     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 435     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 436
 437     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){   //At least common vertex
 438       //    if (( relid1[2]==relid2[2] && coldiff <= 1 )  || ( relid1[3]==relid2[3] &&  rowdiff <= 1 )){ //common side
 439       if((relid1[1] != 0) || CheckTimeGate(d1->GetTime(),d1->GetEnergy(),d2->GetTime(),d2->GetEnergy()))
 440       return 1 ;
 441     }
 442     else {
 443       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
 444         return 2; //  Difference in row numbers is too large to look further
 445     }
 446     return 0 ;
 447
 448   }
 449   else {
 450     if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
 451       return -1 ;
 452     if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
 453       return -1 ;
 454
 455     return 2 ;
 456
 457   }
 458
 459   return 0 ;
 460 }
 461 //____________________________________________________________________________
 462 Bool_t AliPHOSClusterizerv1::CheckTimeGate(Float_t t1, Float_t amp1, Float_t t2, Float_t amp2)const{
 463   //Check if two cells have reasonable time difference
 464   //Note that at low amplitude time is defined up to 1 tick == 100 ns.
 465   if(amp1<fTimeGateLowAmp || amp2<fTimeGateLowAmp){
 466    return (TMath::Abs(t1 - t2 ) < fTimeGateLow) ;
 467   }
 468   else{ //Time should be measured with good accuracy
 469    return (TMath::Abs(t1 - t2 ) < fTimeGateHigh) ;
 470   }
 471
 472 }
 473 //____________________________________________________________________________
 474 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
 475 {
 476   // Tells if (true) or not (false) the digit is in a PHOS-EMC module
 477
 478   Bool_t rv = kFALSE ;
 479
 480   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 481
 482   if(digit->GetId() <= nEMC )   rv = kTRUE;
 483
 484   return rv ;
 485 }
 486
 487 //____________________________________________________________________________
 488 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
 489 {
 490   // Tells if (true) or not (false) the digit is in a PHOS-CPV module
 491
 492   Bool_t rv = kFALSE ;
 493
 494   Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
 495
 496   if(digit->GetId() > nEMC )   rv = kTRUE;
 497
 498   return rv ;
 499 }
 500
 501 //____________________________________________________________________________
 502 void AliPHOSClusterizerv1::WriteRecPoints()
 503 {
 504
 505   // Creates new branches with given title
 506   // fills and writes into TreeR.
 507
 508   Int_t index ;
 509   //Evaluate position, dispersion and other RecPoint properties..
 510   Int_t nEmc = fEMCRecPoints->GetEntriesFast();
 511   Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
 512   TVector3 fakeVtx(0.,0.,0.) ;
 513   for(index = 0; index < nEmc; index++){
 514     AliPHOSEmcRecPoint * rp =
 515       dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
 516     rp->Purify(emcMinE) ;
 517     if(rp->GetMultiplicity()==0){
 518       fEMCRecPoints->RemoveAt(index) ;
 519       delete rp ;
 520       continue;
 521     }
 522
 523     // No vertex is available now, calculate corrections in PID
 524     rp->EvalAll(fDigitsArr) ;
 525     rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
 526     rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
 527     rp->EvalLocal2TrackingCSTransform();
 528   }
 529   fEMCRecPoints->Compress() ;
 530   fEMCRecPoints->Sort() ;
 531   //  fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
 532   for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
 533     dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
 534   }
 535
 536   //For each rec.point set the distance to the nearest bad crystal (BVP)
 537   SetDistancesToBadChannels();
 538
 539   //Now the same for CPV
 540   for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
 541     AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
 542     rp->EvalAll(fDigitsArr) ;
 543     rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
 544     rp->EvalLocal2TrackingCSTransform();
 545   }
 546   fCPVRecPoints->Sort() ;
 547
 548   for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
 549     dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
 550
 551   fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
 552
 553   if(fWrite){ //We write TreeR
 554     fTreeR->Fill();
 555   }
 556 }
 557
 558 //____________________________________________________________________________
 559 void AliPHOSClusterizerv1::MakeClusters()
 560 {
 561   // Steering method to construct the clusters stored in a list of Reconstructed Points
 562   // A cluster is defined as a list of neighbour digits
 563
 564   fNumberOfCpvClusters     = 0 ;
 565   fNumberOfEmcClusters     = 0 ;
 566
 567   //Mark all digits as unused yet
 568   const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;)
 569   Int_t nDigits=fDigitsArr->GetEntriesFast() ;
 570
 571   for(Int_t i=0; i<nDigits; i++){
 572     fDigitsUsed[i]=0 ;
 573   }
 574   Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
 575                      //e.g. first digit in this module, first CPV digit etc.
 576   AliPHOSDigit * digit ;
 577   TArrayI clusterdigitslist(maxNDigits) ;
 578   AliPHOSRecPoint * clu = 0 ;
 579   for(Int_t i=0; i<nDigits; i++){
 580     if(fDigitsUsed[i])
 581       continue ;
 582
 583     digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
 584
 585     clu=0 ;
 586
 587     Int_t index ;
 588
 589     //is this digit so energetic that start cluster?
 590     if (( IsInEmc(digit) &&  Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
 591         ( IsInCpv(digit) &&  Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) {
 592       Int_t iDigitInCluster = 0 ;
 593       if  ( IsInEmc(digit) ) {
 594         // start a new EMC RecPoint
 595         if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 596           fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
 597
 598         fEMCRecPoints->AddAt(new  AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
 599         clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 600         fNumberOfEmcClusters++ ;
 601         clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),CalibrateT(digit->GetTime(),digit->GetId())) ;
 602         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
 603         iDigitInCluster++ ;
 604         fDigitsUsed[i]=kTRUE ;
 605       } else {
 606         // start a new CPV cluster
 607         if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 608           fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
 609
 610         fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
 611         clu =  static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 612         fNumberOfCpvClusters++ ;
 613         clu->AddDigit(*digit,  Calibrate(digit->GetEnergy(),digit->GetId()),0.) ; // no timing information in CPV
 614         clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ;
 615         iDigitInCluster++ ;
 616         fDigitsUsed[i]=kTRUE ;
 617       } // else
 618
 619       //Now scan remaining digits in list to find neigbours of our seed
 620
 621       AliPHOSDigit * digitN ;
 622       index = 0 ;
 623       while (index < iDigitInCluster){ // scan over digits already in cluster
 624         digit =  static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) )  ;
 625         index++ ;
 626         for(Int_t j=iFirst; j<nDigits; j++){
 627           if (iDigitInCluster >= maxNDigits) {
 628             AliError(Form("The number of digits in cluster is more than %d, skip the rest of event",
 629                           maxNDigits));
 630             return;
 631           }
 632           if(fDigitsUsed[j])
 633             continue ;        //look through remaining digits
 634           digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
 635           Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
 636           switch (ineb ) {
 637           case -1:   //too early (e.g. previous module), do not look before j at subsequent passes
 638             iFirst=j ;
 639             break ;
 640           case 0 :   // not a neighbour
 641             break ;
 642           case 1 :   // are neighbours
 643             clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digitN->GetId()),CalibrateT(digitN->GetTime(),digitN->GetId())) ;
 644             clusterdigitslist[iDigitInCluster] = j ;
 645             iDigitInCluster++ ;
 646             fDigitsUsed[j]=kTRUE ;
 647             break ;
 648           case 2 :   // too far from each other
 649             goto endOfLoop;
 650           } // switch
 651
 652         }
 653
 654         endOfLoop: ; //scanned all possible neighbours for this digit
 655
 656       } // loop over cluster
 657     } // energy theshold
 658   }
 659
 660 }
 661
 662 //____________________________________________________________________________
 663 void AliPHOSClusterizerv1::MakeUnfolding()
 664 {
 665   // Unfolds clusters using the shape of an ElectroMagnetic shower
 666   // Performs unfolding of all EMC/CPV clusters
 667
 668   // Unfold first EMC clusters
 669   if(fNumberOfEmcClusters > 0){
 670
 671     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 672
 673     Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
 674     Int_t index ;
 675     for(index = 0 ; index < numberofNotUnfolded ; index++){
 676
 677       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
 678       if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
 679         break ;
 680
 681       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 682       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 683       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 684       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
 685
 686       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 687         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 688
 689         fEMCRecPoints->Remove(emcRecPoint);
 690         fEMCRecPoints->Compress() ;
 691         index-- ;
 692         fNumberOfEmcClusters -- ;
 693         numberofNotUnfolded-- ;
 694       }
 695       else{
 696         emcRecPoint->SetNExMax(1) ; //Only one local maximum
 697       }
 698
 699       delete[] maxAt ;
 700       delete[] maxAtEnergy ;
 701     }
 702   }
 703   // Unfolding of EMC clusters finished
 704
 705
 706   // Unfold now CPV clusters
 707   if(fNumberOfCpvClusters > 0){
 708
 709     Int_t nModulesToUnfold = fGeom->GetNModules() ;
 710
 711     Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
 712     Int_t index ;
 713     for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
 714
 715       AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
 716
 717       if(recPoint->GetPHOSMod()> nModulesToUnfold)
 718         break ;
 719
 720       AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
 721
 722       Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
 723       AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
 724       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 725       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
 726
 727       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 728         UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
 729         fCPVRecPoints->Remove(emcRecPoint);
 730         fCPVRecPoints->Compress() ;
 731         index-- ;
 732         numberofCpvNotUnfolded-- ;
 733         fNumberOfCpvClusters-- ;
 734       }
 735
 736       delete[] maxAt ;
 737       delete[] maxAtEnergy ;
 738     }
 739   }
 740   //Unfolding of Cpv clusters finished
 741
 742 }
 743
 744 //____________________________________________________________________________
 745 Double_t  AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
 746 {
 747   // Shape of the shower (see PHOS TDR)
 748   // If you change this function, change also the gradient evaluation in ChiSquare()
 749
 750   //for the moment we neglect dependence on the incident angle.
 751
 752   Double_t r2    = x*x + z*z ;
 753   Double_t r4    = r2*r2 ;
 754   Double_t r295  = TMath::Power(r2, 2.95/2.) ;
 755   Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
 756   return shape ;
 757 }
 758
 759 //____________________________________________________________________________
 760 void  AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
 761                                           Int_t nMax,
 762                                           AliPHOSDigit ** maxAt,
 763                                           Float_t * maxAtEnergy)
 764 {
 765   // Performs the unfolding of a cluster with nMax overlapping showers
 766
 767   Int_t nPar = 3 * nMax ;
 768   Float_t * fitparameters = new Float_t[nPar] ;
 769
 770   Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
 771
 772   if( !rv ) {
 773     // Fit failed, return and remove cluster
 774     iniEmc->SetNExMax(-1) ;
 775     delete[] fitparameters ;
 776     return ;
 777   }
 778
 779   // create ufolded rec points and fill them with new energy lists
 780   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
 781   // and later correct this number in acordance with actual energy deposition
 782
 783   Int_t nDigits = iniEmc->GetMultiplicity() ;
 784   Float_t * efit = new Float_t[nDigits] ;
 785   Float_t xDigit=0.,zDigit=0. ;
 786   Float_t xpar=0.,zpar=0.,epar=0.  ;
 787   Int_t relid[4] ;
 788   AliPHOSDigit * digit = 0 ;
 789   Int_t * emcDigits = iniEmc->GetDigitsList() ;
 790
 791   TVector3 vIncid ;
 792
 793   Int_t iparam ;
 794   Int_t iDigit ;
 795   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 796     digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
 797     fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 798     fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 799     efit[iDigit] = 0;
 800
 801     iparam = 0 ;
 802     while(iparam < nPar ){
 803       xpar = fitparameters[iparam] ;
 804       zpar = fitparameters[iparam+1] ;
 805       epar = fitparameters[iparam+2] ;
 806       iparam += 3 ;
 807 //      fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
 808 //      efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
 809       efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 810     }
 811   }
 812
 813   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
 814   // so that energy deposited in each cell is distributed betwin new clusters proportionally
 815   // to its contribution to efit
 816
 817   Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
 818   Float_t ratio ;
 819
 820   iparam = 0 ;
 821   while(iparam < nPar ){
 822     xpar = fitparameters[iparam] ;
 823     zpar = fitparameters[iparam+1] ;
 824     epar = fitparameters[iparam+2] ;
 825     iparam += 3 ;
 826 //    fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
 827
 828     AliPHOSEmcRecPoint * emcRP = 0 ;
 829
 830     if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
 831
 832       if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
 833         fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
 834
 835       (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
 836       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
 837       fNumberOfEmcClusters++ ;
 838       emcRP->SetNExMax((Int_t)nPar/3) ;
 839     }
 840     else{//create new entries in fCPVRecPoints
 841       if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
 842         fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
 843
 844       (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
 845       emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
 846       fNumberOfCpvClusters++ ;
 847     }
 848
 849     Float_t eDigit ;
 850     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
 851       digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
 852       fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 853       fGeom->RelPosInModule(relid, xDigit, zDigit) ;
 854 //      ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
 855       ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
 856       eDigit = emcEnergies[iDigit] * ratio ;
 857       emcRP->AddDigit( *digit, eDigit,CalibrateT(digit->GetTime(),digit->GetId()) ) ;
 858     }
 859   }
 860
 861   delete[] fitparameters ;
 862   delete[] efit ;
 863
 864 }
 865
 866 //_____________________________________________________________________________
 867 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
 868 {
 869   // Calculates the Chi square for the cluster unfolding minimization
 870   // Number of parameters, Gradient, Chi squared, parameters, what to do
 871
 872   TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
 873
 874   AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) )  ;
 875   TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
 876   // A bit buggy way to get an access to the geometry
 877   // To be revised!
 878   AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
 879
 880 //  TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ;  //Vertex position
 881
 882   //  AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
 883
 884   Int_t * emcDigits     = emcRP->GetDigitsList() ;
 885
 886   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
 887
 888   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
 889
 890 //  TVector3 vInc ;
 891   fret = 0. ;
 892   Int_t iparam ;
 893
 894   if(iflag == 2)
 895     for(iparam = 0 ; iparam < nPar ; iparam++)
 896       Grad[iparam] = 0 ; // Will evaluate gradient
 897
 898   Double_t efit ;
 899
 900   AliPHOSDigit * digit ;
 901   Int_t iDigit ;
 902
 903   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
 904
 905     digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
 906
 907     Int_t relid[4] ;
 908     Float_t xDigit ;
 909     Float_t zDigit ;
 910
 911     geom->AbsToRelNumbering(digit->GetId(), relid) ;
 912
 913     geom->RelPosInModule(relid, xDigit, zDigit) ;
 914
 915      if(iflag == 2){  // calculate gradient
 916        Int_t iParam = 0 ;
 917        efit = 0 ;
 918        while(iParam < nPar ){
 919          Double_t dx = (xDigit - x[iParam]) ;
 920          iParam++ ;
 921          Double_t dz = (zDigit - x[iParam]) ;
 922          iParam++ ;
 923 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
 924 //         efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
 925          efit += x[iParam] * ShowerShape(dx,dz) ;
 926          iParam++ ;
 927        }
 928        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
 929        iParam = 0 ;
 930        while(iParam < nPar ){
 931          Double_t xpar = x[iParam] ;
 932          Double_t zpar = x[iParam+1] ;
 933          Double_t epar = x[iParam+2] ;
 934 //         fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
 935          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
 936 //         Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
 937          Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 938 //DP: No incident angle dependence in gradient yet!!!!!!
 939          Double_t r4 = dr*dr*dr*dr ;
 940          Double_t r295 = TMath::Power(dr,2.95) ;
 941          Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
 942                                          0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
 943
 944          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
 945          iParam++ ;
 946          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
 947          iParam++ ;
 948          Grad[iParam] += shape ;                                  // Derivative over energy
 949          iParam++ ;
 950        }
 951      }
 952      efit = 0;
 953      iparam = 0 ;
 954
 955      while(iparam < nPar ){
 956        Double_t xpar = x[iparam] ;
 957        Double_t zpar = x[iparam+1] ;
 958        Double_t epar = x[iparam+2] ;
 959        iparam += 3 ;
 960 //       fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
 961 //       efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
 962        efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
 963      }
 964
 965      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
 966      // Here we assume, that sigma = sqrt(E)
 967   }
 968
 969 }
 970
 971 //____________________________________________________________________________
 972 void AliPHOSClusterizerv1::Print(const Option_t *)const
 973 {
 974   // Print clusterizer parameters
 975
 976   TString message ;
 977   TString taskName(GetName()) ;
 978   taskName.ReplaceAll(Version(), "") ;
 979
 980   if( strcmp(GetName(), "") !=0 ) {
 981     // Print parameters
 982     message  = "\n--------------- %s %s -----------\n" ;
 983     message += "Clusterizing digits from the file: %s\n" ;
 984     message += "                           Branch: %s\n" ;
 985     message += "                       EMC Clustering threshold = %f\n" ;
 986     message += "                       EMC Local Maximum cut    = %f\n" ;
 987     message += "                       EMC Logarothmic weight   = %f\n" ;
 988     message += "                       CPV Clustering threshold = %f\n" ;
 989     message += "                       CPV Local Maximum cut    = %f\n" ;
 990     message += "                       CPV Logarothmic weight   = %f\n" ;
 991     if(fToUnfold)
 992       message += " Unfolding on\n" ;
 993     else
 994       message += " Unfolding off\n" ;
 995
 996     message += "------------------------------------------------------------------" ;
 997   }
 998   else
 999     message = " AliPHOSClusterizerv1 not initialized " ;
1000
1001   AliInfo(Form("%s, %s %s %s %s %f %f %f %f %f %f", message.Data(),
1002        taskName.Data(),
1003        GetTitle(),
1004        taskName.Data(),
1005        GetName(),
1006        fEmcClusteringThreshold,
1007        fEmcLocMaxCut,
1008        fW0,
1009        fCpvClusteringThreshold,
1010        fCpvLocMaxCut,
1011        fW0CPV )) ;
1012 }
1013 //____________________________________________________________________________
1014 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1015 {
1016   // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1017
1018   AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
1019                fEMCRecPoints->GetEntriesFast(),
1020                fCPVRecPoints->GetEntriesFast() ))  ;
1021
1022   if(strstr(option,"all")) {
1023     printf("\n EMC clusters \n") ;
1024     printf("Index    Ene(MeV) Multi Module     X    Y   Z    Lambdas_1  Lambda_2  # of prim  Primaries list\n") ;
1025     Int_t index ;
1026     for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1027       AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1028       TVector3  locpos;
1029       rp->GetLocalPosition(locpos);
1030       Float_t lambda[2];
1031       rp->GetElipsAxis(lambda);
1032       Int_t * primaries;
1033       Int_t nprimaries;
1034       primaries = rp->GetPrimaries(nprimaries);
1035       printf("\n%6d  %8.2f  %3d     %2d     %4.1f   %4.1f   %4.1f   %4f  %4f    %2d     : ",
1036               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1037               locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1038
1039       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1040         printf("%d ", primaries[iprimary] ) ;
1041       }
1042       printf("\n") ;
1043     }
1044
1045     //Now plot CPV recPoints
1046     printf("\n CPV clusters \n") ;
1047     printf("Index    Ene(MeV) Module     X     Y    Z  \n") ;
1048     for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1049       AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1050
1051       TVector3  locpos;
1052       rp->GetLocalPosition(locpos);
1053
1054       printf("\n%6d  %8.2f  %2d     %4.1f    %4.1f %4.1f \n",
1055              rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1056              locpos.X(), locpos.Y(), locpos.Z()) ;
1057     }
1058   }
1059 }
1060
1061
1062 //____________________________________________________________________________
1063 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1064 {
1065   //For each EMC rec. point set the distance to the nearest bad crystal.
1066   //Author: Boris Polichtchouk
1067
1068   if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1069
1070   Int_t badIds[8000];
1071   memset(badIds,0,8000*sizeof(Int_t));
1072   fgCalibData->EmcBadChannelIds(badIds);
1073
1074   AliPHOSEmcRecPoint* rp;
1075
1076   TMatrixF gmat;
1077   TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1078   TVector3 gposBadChannel; // global position of bad crystal
1079   TVector3 dR;
1080
1081   Float_t dist,minDist;
1082   Int_t relid[4]={0,0,0,0} ;
1083   TVector3 lpos ;
1084   for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1085     rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1086     //evaluate distance to border
1087     relid[0]=rp->GetPHOSMod() ;
1088     relid[2]=1 ;
1089     relid[3]=1 ;
1090     Float_t xcorner,zcorner;
1091     fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1092     rp->GetLocalPosition(lpos) ;
1093     minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1094     for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1095       fGeom->AbsToRelNumbering(badIds[iBad],relid)  ;
1096       if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1097         continue ;                   //bad channels can be in the module which does not exist in simulations.
1098       rp->GetGlobalPosition(gposRecPoint,gmat);
1099       fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1100       AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1101                       gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1102                       gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1103       dR = gposBadChannel-gposRecPoint;
1104       dist = dR.Mag();
1105       if(dist<minDist) minDist = dist;
1106     }
1107
1108     rp->SetDistanceToBadCrystal(minDist);
1109   }
1110
1111 }
1112 //==================================================================================
1113 Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId) const{
1114   // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
1115
1116   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1117
1118   //Determine rel.position of the cell absolute ID
1119   Int_t relId[4];
1120   geom->AbsToRelNumbering(absId,relId);
1121   Int_t module=relId[0];
1122   Int_t row   =relId[2];
1123   Int_t column=relId[3];
1124   if(relId[1]){ //CPV
1125     Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row);
1126     return amp*calibration ;
1127   }
1128   else{ //EMC
1129     Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row);
1130     return amp*calibration ;
1131   }
1132 }
1133 //==================================================================================
1134 Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId)const{
1135   // Calibrate time in EMC digit
1136
1137   const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1138
1139   //Determine rel.position of the cell absolute ID
1140   Int_t relId[4];
1141   geom->AbsToRelNumbering(absId,relId);
1142   Int_t module=relId[0];
1143   Int_t row   =relId[2];
1144   Int_t column=relId[3];
1145   if(relId[1]){ //CPV
1146     return 0. ;
1147   }
1148   else{ //EMC
1149     time += fgCalibData->GetTimeShiftEmc(module,column,row);
1150     return time ;
1151   }
1152 }
1153