Fixed memory leaks for #86360: High memory consumption in 2.76TeV p+p RAW reco jobs
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALReconstructor.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 //_________________________________________________________________________
19 //--
20 //-- Yves Schutz (SUBATECH) 
21 // Reconstruction class. Redesigned from the old AliReconstructionner class and 
22 // derived from STEER/AliReconstructor. 
23 // 
24
25
26 // --- ROOT system ---
27 #include <TClonesArray.h>
28 #include "TGeoManager.h"
29 #include "TGeoMatrix.h"
30
31 // --- Standard library ---
32
33 // --- AliRoot header files ---
34 #include "AliEMCALReconstructor.h"
35
36 #include "AliCodeTimer.h"
37 #include "AliCaloCalibPedestal.h"
38 #include "AliEMCALCalibData.h"
39 #include "AliESDEvent.h"
40 #include "AliESDCaloCluster.h"
41 #include "AliESDCaloCells.h"
42 #include "AliESDtrack.h"
43 #include "AliEMCALLoader.h"
44 #include "AliEMCALRawUtils.h"
45 #include "AliEMCALDigit.h"
46 #include "AliEMCALClusterizerv1.h"
47 #include "AliEMCALClusterizerv2.h"
48 #include "AliEMCALClusterizerNxN.h"
49 #include "AliEMCALRecPoint.h"
50 #include "AliEMCALPID.h"
51 #include "AliRawReader.h"
52 #include "AliCDBEntry.h"
53 #include "AliCDBManager.h"
54 #include "AliEMCALGeometry.h"
55 #include "AliEMCAL.h"
56 #include "AliESDVZERO.h"
57 #include "AliCDBManager.h"
58 #include "AliRunLoader.h"
59 #include "AliRun.h"
60 #include "AliEMCALTriggerData.h"
61 #include "AliEMCALTriggerElectronics.h"
62 #include "AliEMCALTriggerDCSConfigDB.h"
63 #include "AliEMCALTriggerDCSConfig.h"
64 #include "AliEMCALTriggerData.h"
65 #include "AliEMCALTriggerRawDigit.h"
66 #include "AliEMCALTriggerPatch.h"
67 #include "AliEMCALTriggerTypes.h"
68
69 ClassImp(AliEMCALReconstructor) 
70   
71 const AliEMCALRecParam*     AliEMCALReconstructor::fgkRecParam        = 0;   // EMCAL rec. parameters
72 AliEMCALRawUtils*           AliEMCALReconstructor::fgRawUtils         = 0;   // EMCAL raw utilities class
73 AliEMCALClusterizer*        AliEMCALReconstructor::fgClusterizer      = 0;   // EMCAL clusterizer class
74 TClonesArray*               AliEMCALReconstructor::fgDigitsArr        = 0;   // list of digits, to be used multiple times
75 TObjArray*                  AliEMCALReconstructor::fgClustersArr      = 0;   // list of clusters, to be used multiple times
76 TClonesArray*               AliEMCALReconstructor::fgTriggerDigits    = 0;   // list of trigger digits, to be used multiple times
77 AliEMCALTriggerElectronics* AliEMCALReconstructor::fgTriggerProcessor = 0x0;
78 //____________________________________________________________________________
79 AliEMCALReconstructor::AliEMCALReconstructor() 
80   : fGeom(0),fCalibData(0),fPedestalData(0),fTriggerData(0x0), fMatches(0x0)
81 {
82   // ctor
83
84   // AliDebug(2, "Mark.");  
85
86   fgRawUtils = new AliEMCALRawUtils;
87   
88   //To make sure we match with the geometry in a simulation file,
89   //let's try to get it first.  If not, take the default geometry
90   AliRunLoader *rl = AliRunLoader::Instance();
91   if (rl->GetAliRun()){
92     AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
93     if(emcal) fGeom = emcal->GetGeometry();
94   }
95   
96   if(!fGeom) {
97     AliInfo(Form("Using default geometry in reconstruction"));
98     fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
99   }
100   
101   //Get calibration parameters  
102   if(!fCalibData)
103     {
104       AliCDBEntry *entry = (AliCDBEntry*) 
105         AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
106       if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
107     }
108   
109   if(!fCalibData)
110     AliFatal("Calibration parameters not found in CDB!");
111   
112   //Get calibration parameters  
113   if(!fPedestalData)
114     {
115       AliCDBEntry *entry = (AliCDBEntry*) 
116         AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
117       if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
118     }
119   
120   if(!fPedestalData)
121     AliFatal("Dead map not found in CDB!");
122   
123   if(!fGeom) AliFatal(Form("Could not get geometry!"));
124   
125   AliEMCALTriggerDCSConfigDB* dcsConfigDB = AliEMCALTriggerDCSConfigDB::Instance();
126   
127   const AliEMCALTriggerDCSConfig* dcsConfig = dcsConfigDB->GetTriggerDCSConfig();
128   
129   if (!dcsConfig) AliFatal("No Trigger DCS Configuration from OCDB!");
130   fgTriggerProcessor = new AliEMCALTriggerElectronics( dcsConfig );
131   
132   fTriggerData = new AliEMCALTriggerData();
133   
134   //Init temporary list of digits
135   fgDigitsArr     = new TClonesArray("AliEMCALDigit",1000);
136   fgClustersArr   = new TObjArray(1000);
137   fgTriggerDigits = new TClonesArray("AliEMCALTriggerRawDigit",1000);   
138
139   //Track matching
140   fMatches = new TList();
141   fMatches->SetOwner(kTRUE);
142
143
144 //____________________________________________________________________________
145 AliEMCALReconstructor::~AliEMCALReconstructor()
146 {
147   // dtor
148
149   //AliDebug(2, "Mark.");
150
151   if(fGeom)              delete fGeom;
152   
153   //No need to delete, recovered from OCDB
154   //if(fCalibData)         delete fCalibData;
155   //if(fPedestalData)      delete fPedestalData;
156   
157   if(fgDigitsArr){
158     fgDigitsArr->Clear("C");
159     delete fgDigitsArr; 
160   }
161   
162   if(fgClustersArr){
163     fgClustersArr->Clear();
164     delete fgClustersArr; 
165   }
166   
167   if(fgTriggerDigits){
168     fgTriggerDigits->Clear();
169     delete fgTriggerDigits; 
170   }
171   
172   if(fgRawUtils)         delete fgRawUtils;
173   if(fgClusterizer)      delete fgClusterizer;
174   if(fgTriggerProcessor) delete fgTriggerProcessor;
175   
176   if(fMatches) { fMatches->Delete(); delete fMatches; fMatches = 0;}
177   
178   AliCodeTimer::Instance()->Print();
179
180
181 //____________________________________________________________________________                                  
182 void AliEMCALReconstructor::InitClusterizer() const
183 {
184   //Init the clusterizer with geometry and calibration pointers, avoid doing it twice.                          
185   Int_t clusterizerType = -1;
186   Int_t eventType = -1;
187   if(GetRecParam()) {
188     clusterizerType = GetRecParam()->GetClusterizerFlag();
189     eventType       = GetRecParam()->GetEventSpecie();
190   }
191   else{
192     AliCDBEntry *entry = (AliCDBEntry*)
193       AliCDBManager::Instance()->Get("EMCAL/Calib/RecoParam");
194     //Get The reco param for the default event specie                                                           
195     if (entry) {
196       AliEMCALRecParam *recParam  = (AliEMCALRecParam*)((TObjArray *) entry->GetObject())->At(0);
197       if(recParam) clusterizerType = recParam->GetClusterizerFlag(); 
198     }
199   }
200   
201   //Check if clusterizer previously set corresponds to what is needed for this event type                       
202   if(fgClusterizer){
203     if(eventType!=AliRecoParam::kCalib){
204       //printf("ReCreate clusterizer? Clusterizer set <%d>, Clusterizer in use <%s>\n",
205       //     clusterizerType, fgClusterizer->Version());
206       
207       if     (clusterizerType == AliEMCALRecParam::kClusterizerv1  && !strcmp(fgClusterizer->Version(),"clu-v1"))  return;
208       
209       else if(clusterizerType == AliEMCALRecParam::kClusterizerNxN && !strcmp(fgClusterizer->Version(),"clu-NxN")) return;
210       
211       else if(clusterizerType == AliEMCALRecParam::kClusterizerv2  && !strcmp(fgClusterizer->Version(),"clu-v2"))  return;
212       
213       //Need to create new clusterizer, the one set previously is not the correct one     
214       delete fgClusterizer;
215     }
216     else return;
217   }
218   
219   if      (clusterizerType  == AliEMCALRecParam::kClusterizerv1)
220     {
221       fgClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData,fPedestalData);
222     }
223   else if (clusterizerType  == AliEMCALRecParam::kClusterizerNxN)
224     {
225       fgClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData,fPedestalData);
226     }
227   else if (clusterizerType  == AliEMCALRecParam::kClusterizerv2)
228   {
229     fgClusterizer = new AliEMCALClusterizerv2   (fGeom, fCalibData,fPedestalData);
230   }
231   else 
232   {
233     AliFatal(Form("Unknown clusterizer %d ", clusterizerType));
234   }
235 }
236
237 //____________________________________________________________________________
238 void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
239 {
240   // method called by AliReconstruction; 
241   // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
242   // segment maker needs access to the AliESD object to retrieve the tracks reconstructed by 
243   // the global tracking.
244   // Works on the current event.
245   
246   AliCodeTimerAuto("",0)
247     
248   //Get input digits and put them in fgDigitsArr, clear the list before 
249   ReadDigitsArrayFromTree(digitsTree);
250   
251   InitClusterizer();
252   
253   fgClusterizer->InitParameters();
254   fgClusterizer->SetOutput(clustersTree);
255   
256   //Skip clusterization of LED events
257   if (GetRecParam()->GetEventSpecie()!=AliRecoParam::kCalib){
258     
259     if(fgDigitsArr && fgDigitsArr->GetEntries()) {
260       
261       fgClusterizer->SetInput(digitsTree);
262       
263       //fgClusterizer->Digits2Clusters("deb all") ; //For debugging
264       fgClusterizer->Digits2Clusters("");
265       
266       fgClusterizer->Clear();
267       
268     }//digits array exists and has somethind
269   }//not a LED event
270   
271   clustersTree->Fill(); 
272
273   // Deleting the recpoints at the end of the reconstruction call
274   fgClusterizer->DeleteRecPoints();
275 }
276
277 //____________________________________________________________________________
278 void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
279   
280 {
281   // Conversion from raw data to
282   // EMCAL digits.
283   // Works on a single-event basis
284   
285   rawReader->Reset() ; 
286   
287   fTriggerData->SetMode(1);     
288   
289   if(fgDigitsArr) fgDigitsArr->Clear("C");
290   
291   TClonesArray *digitsTrg = new TClonesArray("AliEMCALTriggerRawDigit", 32 * 96);
292   
293   Int_t bufsize = 32000;
294   digitsTree->Branch("EMCAL", &fgDigitsArr, bufsize);
295   digitsTree->Branch("EMTRG", &digitsTrg, bufsize);
296   
297   //Skip calibration events do the rest
298   Bool_t doFit = kTRUE;
299   if ( !(GetRecParam()->FitLEDEvents()) && GetRecParam()->GetEventSpecie()==AliRecoParam::kCalib) doFit = kFALSE;
300   if (doFit){
301     //must be done here because, in constructor, option is not yet known
302     fgRawUtils->SetOption(GetOption());
303     
304     //  fgRawUtils->SetRawFormatHighLowGainFactor(GetRecParam()->GetHighLowGainFactor());
305     
306     //   fgRawUtils->SetRawFormatOrder(GetRecParam()->GetOrderParameter());
307     //    fgRawUtils->SetRawFormatTau(GetRecParam()->GetTau());
308     fgRawUtils->SetNoiseThreshold(GetRecParam()->GetNoiseThreshold());
309     fgRawUtils->SetNPedSamples(GetRecParam()->GetNPedSamples());
310     fgRawUtils->SetRemoveBadChannels(GetRecParam()->GetRemoveBadChannels());
311     if (!fgRawUtils->GetFittingAlgorithm()) fgRawUtils->SetFittingAlgorithm(GetRecParam()->GetFittingAlgorithm());
312     fgRawUtils->SetFALTROUsage(GetRecParam()->UseFALTRO());
313     
314     //fgRawUtils->SetTimeMin(GetRecParam()->GetTimeMin());
315     //fgRawUtils->SetTimeMax(GetRecParam()->GetTimeMax());
316     
317     //  fgRawUtils->SetTimeMin(-99999 );
318     //  fgRawUtils->SetTimeMax( 99999 );
319     
320     fgRawUtils->Raw2Digits(rawReader,fgDigitsArr,fPedestalData,digitsTrg,fTriggerData);
321     
322   }//skip calibration event
323   else{
324     AliDebug(1," Calibration Event, skip!");
325   }
326   
327   digitsTree->Fill();
328   digitsTrg->Delete();
329   delete digitsTrg;
330   
331 }
332
333
334 //____________________________________________________________________________
335 void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, 
336                                     AliESDEvent* esd) const
337 {
338   // Called by AliReconstruct after Reconstruct() and global tracking and vertexing 
339   // and V0 
340   // Works on the current event
341   // printf(" ## AliEMCALReconstructor::FillESD() is started ### \n ");
342   //return;
343   
344   //########################################
345   // Trigger
346   //########################################
347   
348   static int saveOnce = 0;
349   
350   Int_t v0M[2] = {0, 0};
351   
352   AliESDVZERO* esdV0 = esd->GetVZEROData();
353   
354   if (esdV0) 
355     {
356       for (Int_t i = 0; i < 32; i++)
357         {
358           v0M[0] += (Int_t)esdV0->GetAdcV0C(i);
359           v0M[1] += (Int_t)esdV0->GetAdcV0A(i);
360         }
361     }
362   else
363     {
364       AliWarning("Cannot retrieve V0 ESD! Run w/ null V0 charges");
365     }
366   
367   if (fgTriggerDigits) fgTriggerDigits->Clear();
368   
369   TBranch *branchtrg = digitsTree->GetBranch("EMTRG");
370   
371   if (!branchtrg) 
372     { 
373       AliError("Can't get the branch with the EMCAL trigger digits!");
374       return;
375     }
376   
377   branchtrg->SetAddress(&fgTriggerDigits);
378   branchtrg->GetEntry(0);
379   
380   // Note: fgTriggerProcessor reset done at the end of this method
381   fgTriggerProcessor->Digits2Trigger(fgTriggerDigits, v0M, fTriggerData);
382   
383   // Fill ESD
384   AliESDCaloTrigger* trgESD = esd->GetCaloTrigger("EMCAL");
385   
386   if (trgESD)
387     {
388       trgESD->Allocate(fgTriggerDigits->GetEntriesFast());
389       
390       for (Int_t i = 0; i < fgTriggerDigits->GetEntriesFast(); i++)
391         {         
392           AliEMCALTriggerRawDigit* rdig = (AliEMCALTriggerRawDigit*)fgTriggerDigits->At(i);
393           
394           Int_t px, py;
395           if (fGeom->GetPositionInEMCALFromAbsFastORIndex(rdig->GetId(), px, py))
396             {
397               Int_t a = -1, t = -1, times[10]; 
398               
399               rdig->GetMaximum(a, t);
400               rdig->GetL0Times(times);
401               
402               trgESD->Add(px, py, a, t, times, rdig->GetNL0Times(), rdig->GetL1TimeSum(), rdig->GetTriggerBits());
403             }
404         }
405       
406       trgESD->SetL1Threshold(0, fTriggerData->GetL1GammaThreshold());
407       
408       trgESD->SetL1Threshold(1, fTriggerData->GetL1JetThreshold()  );
409       
410       Int_t v0[2];
411       fTriggerData->GetL1V0(v0);
412       
413       trgESD->SetL1V0(v0);      
414       trgESD->SetL1FrameMask(fTriggerData->GetL1FrameMask());            
415       
416       if (!saveOnce && fTriggerData->GetL1DataDecoded()) 
417         {
418           int type[8] = {0};
419           fTriggerData->GetL1TriggerType(type);
420           
421           esd->SetCaloTriggerType(type);
422           
423           saveOnce = 1;
424         }
425     }
426   
427   // Resetting
428   fTriggerData->Reset();
429   
430   //########################################
431   //##############Fill CaloCells###############
432   //########################################
433   
434   //Get input digits and put them in fgDigitsArr, clear the list before 
435   ReadDigitsArrayFromTree(digitsTree);
436   
437   Int_t nDigits = fgDigitsArr->GetEntries(), idignew = 0 ;
438   AliDebug(1,Form("%d digits",nDigits));
439   AliESDCaloCells &emcCells = *(esd->GetEMCALCells());
440   emcCells.CreateContainer(nDigits);
441   emcCells.SetType(AliVCaloCells::kEMCALCell);
442   Float_t energy = 0;
443   Float_t time   = 0;
444   for (Int_t idig = 0 ; idig < nDigits ; idig++) {
445     const AliEMCALDigit * dig = (const AliEMCALDigit*)fgDigitsArr->At(idig);
446     time   = dig->GetTime();      // Time already calibrated in clusterizer
447     energy = dig->GetAmplitude(); // energy calibrated in clusterizer
448     if(energy > 0 ){
449       fgClusterizer->Calibrate(energy,time,dig->GetId()); //Digits already calibrated in clusterizers
450       if(energy > 0){ //Digits tagged as bad (dead, hot, not alive) are set to 0 in calibrate, remove them
451         emcCells.SetCell(idignew,dig->GetId(),energy, time);   
452         idignew++;
453       }
454     }
455   }
456   emcCells.SetNumberOfCells(idignew);
457   emcCells.Sort();
458   
459   //------------------------------------------------------------
460   //-----------------CLUSTERS-----------------------------
461   //------------------------------------------------------------
462   clustersTree->SetBranchStatus("*",0); //disable all branches
463   clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
464   if(fgClustersArr) fgClustersArr->Clear();
465   TBranch *branch = clustersTree->GetBranch("EMCALECARP");
466   branch->SetAddress(&fgClustersArr);
467   branch->GetEntry(0);
468   //clustersTree->GetEvent(0);
469   
470   Int_t nClusters = fgClustersArr->GetEntries(),  nClustersNew=0;
471   AliDebug(1,Form("%d clusters",nClusters));
472
473   
474   //########################################
475   //##############Fill CaloClusters#############
476   //########################################
477   for (Int_t iClust = 0 ; iClust < nClusters ; iClust++) {
478     const AliEMCALRecPoint * clust = (const AliEMCALRecPoint*)fgClustersArr->At(iClust);
479     if(!clust) continue;
480     //if(clust->GetClusterType()== AliVCluster::kEMCALClusterv1) nRP++; else nPC++;
481     // clust->Print(); //For debugging
482     // Get information from EMCAL reconstruction points
483     Float_t xyz[3];
484     TVector3 gpos;
485     clust->GetGlobalPosition(gpos);
486     for (Int_t ixyz=0; ixyz<3; ixyz++)
487       xyz[ixyz] = gpos[ixyz];
488     Float_t elipAxis[2];
489     clust->GetElipsAxis(elipAxis);
490     //Create digits lists
491     Int_t cellMult = clust->GetMultiplicity();
492     //TArrayS digiList(digitMult);
493     Float_t *amplFloat = clust->GetEnergiesList();
494     Int_t   *digitInts = clust->GetAbsId();
495     TArrayS absIdList(cellMult);
496     TArrayD fracList(cellMult);
497     
498     Int_t newCellMult = 0;
499     for (Int_t iCell=0; iCell<cellMult; iCell++) {
500       if (amplFloat[iCell] > 0) {
501         absIdList[newCellMult] = (UShort_t)(digitInts[iCell]);
502         //Calculate Fraction
503         if(emcCells.GetCellAmplitude(digitInts[iCell])>0 && GetRecParam()->GetUnfold()){
504           fracList[newCellMult] = amplFloat[iCell]/(emcCells.GetCellAmplitude(digitInts[iCell]));//get cell calibration value 
505           
506         }
507         else{
508           fracList[newCellMult] = 0; 
509         }
510         newCellMult++;
511       }
512     }
513     
514     absIdList.Set(newCellMult);
515     fracList.Set(newCellMult);
516     
517     if(newCellMult > 0) { // accept cluster if it has some digit
518       nClustersNew++;
519       //Primaries
520       Int_t  parentMult  = 0;
521       Int_t *parentList =  clust->GetParents(parentMult);
522       // fills the ESDCaloCluster
523       AliESDCaloCluster * ec = new AliESDCaloCluster() ;
524       ec->SetType(AliVCluster::kEMCALClusterv1);
525       ec->SetPosition(xyz);
526       ec->SetE(clust->GetEnergy());
527       
528       //Distance to the nearest bad crystal
529       ec->SetDistanceToBadChannel(clust->GetDistanceToBadTower()); 
530       
531       ec->SetNCells(newCellMult);
532       //Change type of list from short to ushort
533       UShort_t *newAbsIdList  = new UShort_t[newCellMult];
534       Double_t *newFracList   = new Double_t[newCellMult];
535       for(Int_t i = 0; i < newCellMult ; i++) {
536         newAbsIdList[i]=absIdList[i];
537         newFracList[i] =fracList[i];
538       }
539       ec->SetCellsAbsId(newAbsIdList);
540       ec->SetCellsAmplitudeFraction(newFracList);
541       ec->SetDispersion(clust->GetDispersion());
542       ec->SetChi2(-1); //not yet implemented
543       ec->SetM02(elipAxis[0]*elipAxis[0]) ;
544       ec->SetM20(elipAxis[1]*elipAxis[1]) ;
545       ec->SetTOF(clust->GetTime()) ; //time-of-fligh
546       ec->SetNExMax(clust->GetNExMax());          //number of local maxima
547   
548       
549       TArrayI arrayParents(parentMult,parentList);
550       ec->AddLabels(arrayParents);
551       //
552       //Track matching
553       //
554       fMatches->Clear();
555       Int_t nTracks = esd->GetNumberOfTracks();
556       for (Int_t itrack = 0; itrack < nTracks; itrack++)
557         {
558           AliESDtrack * track = esd->GetTrack(itrack) ; // retrieve track
559           if(track->GetEMCALcluster()==iClust)
560             {
561               Double_t dEta=-999, dPhi=-999;
562               Bool_t isMatch =  CalculateResidual(track, ec, dEta, dPhi);
563               if(!isMatch) 
564                 {
565                   // AliDebug(10, "Not good");
566                   continue;
567                 }
568               AliEMCALMatch *match = new AliEMCALMatch();
569               match->SetIndexT(itrack);
570               match->SetDistance(TMath::Sqrt(dEta*dEta+dPhi*dPhi));
571               match->SetdEta(dEta);
572               match->SetdPhi(dPhi);
573               fMatches->Add(match);
574             }
575         } 
576       fMatches->Sort(kSortAscending); //Sort matched tracks from closest to furthest
577       Int_t nMatch = fMatches->GetEntries();
578       TArrayI arrayTrackMatched(nMatch);
579       for(Int_t imatch=0; imatch<nMatch; imatch++)
580         {
581           AliEMCALMatch *match = (AliEMCALMatch*)fMatches->At(imatch);
582           arrayTrackMatched[imatch] = match->GetIndexT();
583           if(imatch==0)
584             {
585               ec->SetTrackDistance(match->GetdPhi(), match->GetdEta());
586             }
587         }
588       ec->AddTracksMatched(arrayTrackMatched);
589     
590       //add the cluster to the esd object
591       esd->AddCaloCluster(ec);
592
593       delete ec;
594       delete [] newAbsIdList ;
595       delete [] newFracList ;
596     }
597   } // cycle on clusters
598
599   //
600   //Reset the index of matched cluster for tracks
601   //to the one in CaloCluster array
602   Int_t ncls = esd->GetNumberOfCaloClusters();
603   for(Int_t icl=0; icl<ncls; icl++)
604     {
605       AliESDCaloCluster *cluster = esd->GetCaloCluster(icl);
606       if(!cluster || !cluster->IsEMCAL()) continue;
607       TArrayI *trackIndex = cluster->GetTracksMatched();
608       for(Int_t itr=0; itr<trackIndex->GetSize(); itr++)
609         {
610           AliESDtrack *track = esd->GetTrack(trackIndex->At(itr));
611           track->SetEMCALcluster(cluster->GetID());
612         }
613     }
614   
615   
616   //Fill ESDCaloCluster with PID weights
617   AliEMCALPID *pid = new AliEMCALPID;
618   //pid->SetPrintInfo(kTRUE);
619   pid->SetReconstructor(kTRUE);
620   pid->RunPID(esd);
621   delete pid;
622   
623   //Store EMCAL misalignment matrixes
624   FillMisalMatrixes(esd) ;
625   
626 }
627
628 //==================================================================================
629 void AliEMCALReconstructor::FillMisalMatrixes(AliESDEvent* esd)const{
630   //Store EMCAL matrixes in ESD Header
631   
632   //Check, if matrixes was already stored
633   for(Int_t sm = 0 ; sm < fGeom->GetNumberOfSuperModules(); sm++){
634     if(esd->GetEMCALMatrix(sm)!=0)
635       return ;
636   }
637   
638   //Create and store matrixes
639   if(!gGeoManager){
640     AliError("Can not store misal. matrixes: no gGeoManager! \n") ;
641     return ;
642   }
643   //Note, that owner of copied marixes will be header
644   const Int_t bufsize = 255;
645   char path[bufsize] ;
646   TGeoHMatrix * m = 0x0;
647   for(Int_t sm = 0; sm < fGeom->GetNumberOfSuperModules(); sm++){
648     snprintf(path,bufsize,"/ALIC_1/XEN1_1/SMOD_%d",sm+1) ; //In Geometry modules numbered 1,2,.,5
649     if(sm >= 10) snprintf(path,bufsize,"/ALIC_1/XEN1_1/SM10_%d",sm-10+1) ;
650     
651     if (gGeoManager->CheckPath(path)){
652       gGeoManager->cd(path);
653       m = gGeoManager->GetCurrentMatrix() ;
654       //                        printf("================================================= \n");
655       //                        printf("AliEMCALReconstructor::FixMisalMatrixes(), sm %d, \n",sm);
656       //                        m->Print("");
657       esd->SetEMCALMatrix(new TGeoHMatrix(*m),sm) ;
658       //                        printf("================================================= \n");
659     }
660     else{
661       esd->SetEMCALMatrix(NULL,sm) ;
662     }
663   }
664 }
665
666 //__________________________________________________________________________
667 void AliEMCALReconstructor::ReadDigitsArrayFromTree(TTree *digitsTree) const
668 {
669   // Read the digits from the input tree
670   // See AliEMCALClusterizer::SetInput(TTree *digitsTree);    
671   
672   // Clear previous digits in the list
673   if(fgDigitsArr){ 
674     fgDigitsArr->Clear("C");
675   }
676   else{
677     // It should not happen, but just in case ...
678     fgDigitsArr = new TClonesArray("AliEMCALDigit",100); 
679   }
680   
681   // Read the digits from the input tree
682   TBranch *branch = digitsTree->GetBranch("EMCAL");
683   if (!branch) { 
684     AliError("can't get the branch with the EMCAL digits !");
685     return;
686   }  
687   
688   branch->SetAddress(&fgDigitsArr);
689   branch->GetEntry(0);
690 }
691
692 //==================================================================================
693 Bool_t AliEMCALReconstructor::CalculateResidual(AliESDtrack *track, AliESDCaloCluster *cluster, Double_t &dEta, Double_t &dPhi)const
694 {
695   //
696   // calculate the residual between track and cluster
697   //
698
699   // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
700   // Otherwise use the TPCInner point
701   AliExternalTrackParam *trkParam = 0;
702   const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
703   if(friendTrack && friendTrack->GetTPCOut())
704     trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
705   else
706     trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
707   if(!trkParam) return kFALSE;
708
709   //Perform extrapolation
710   Double_t trkPos[3];
711   Float_t  clsPos[3];
712
713   AliExternalTrackParam trkParamTmp (*trkParam);
714   cluster->GetPosition(clsPos);
715   TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
716   Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
717   //Rotate to proper local coordinate
718   vec.RotateZ(-alpha); 
719   trkParamTmp.Rotate(alpha); 
720   //extrapolation is done here
721   if(!AliTrackerBase::PropagateTrackToBxByBz(&trkParamTmp, vec.X(), track->GetMass(), GetRecParam()->GetExtrapolateStep(), kFALSE, 0.8, -1)) 
722     return kFALSE; 
723
724   //Calculate the residuals
725   trkParamTmp.GetXYZ(trkPos); 
726    
727   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
728   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
729       
730   Double_t clsPhi = clsPosVec.Phi();
731   if(clsPhi<0) clsPhi+=2*TMath::Pi();
732   Double_t trkPhi = trkPosVec.Phi();
733   if(trkPhi<0) trkPhi+=2*TMath::Pi();
734
735   dPhi = clsPhi-trkPhi;
736   dEta = clsPosVec.Eta()-trkPosVec.Eta();
737
738   return kTRUE;
739 }
740
741 //
742 //==================================================================================
743 //
744 AliEMCALReconstructor::AliEMCALMatch::AliEMCALMatch() 
745   : TObject(),  
746     fIndexT(-1), 
747     fDistance(-999.),
748     fdEta(-999.),
749     fdPhi(-999.)
750 {
751   //default constructor
752
753 }
754
755 //
756 //==================================================================================
757 //
758 AliEMCALReconstructor::AliEMCALMatch::AliEMCALMatch(const AliEMCALMatch& copy)
759   : TObject(),
760     fIndexT(copy.fIndexT),
761     fDistance(copy.fDistance),
762     fdEta(copy.fdEta),
763     fdPhi(copy.fdPhi)
764 {
765   //copy ctor
766 }
767
768 //
769 //==================================================================================
770 //
771 Int_t AliEMCALReconstructor::AliEMCALMatch::Compare(const TObject *obj) const 
772 {
773   //
774   // Compare wrt the residual
775   //
776         
777   AliEMCALReconstructor::AliEMCALMatch *that = (AliEMCALReconstructor::AliEMCALMatch*)obj;
778         
779   Double_t thisDist = fDistance;//fDistance;
780   Double_t thatDist = that->fDistance;//that->GetDistance();
781         
782   if (thisDist > thatDist) return 1;
783   else if (thisDist < thatDist) return -1;
784   return 0;
785 }