modification to read Altro mappings into AliEMCALRecParam and pass to AliEMCALRawUtil...
[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 // --- ROOT system ---
25
26 // --- Standard library ---
27
28 // --- AliRoot header files ---
29 #include "AliEMCALReconstructor.h"
30
31 #include "AliESDEvent.h"
32 #include "AliESDCaloCluster.h"
33 #include "AliESDCaloCells.h"
34 #include "AliESDtrack.h"
35 #include "AliEMCALLoader.h"
36 #include "AliEMCALRawUtils.h"
37 #include "AliEMCALDigit.h"
38 #include "AliEMCALClusterizerv1.h"
39 #include "AliEMCALRecPoint.h"
40 #include "AliEMCALPID.h"
41 #include "AliEMCALTrigger.h"
42 #include "AliRawReader.h"
43 #include "AliCDBEntry.h"
44 #include "AliCDBManager.h"
45 #include "AliEMCALRecParam.h"
46 // to be removed - it is here just because of geom
47 #include "AliRun.h"
48 #include "AliRunLoader.h"
49
50 ClassImp(AliEMCALReconstructor)
51
52 AliEMCALRecParam* AliEMCALReconstructor::fgkRecParam = 0;  // EMCAL rec. parameters
53
54 //____________________________________________________________________________
55 AliEMCALReconstructor::AliEMCALReconstructor() 
56   : fDebug(kFALSE) 
57 {
58   // ctor
59
60
61
62 //____________________________________________________________________________
63 AliEMCALReconstructor::AliEMCALReconstructor(const AliEMCALReconstructor & rec)
64   : AliReconstructor(rec),
65     fDebug(rec.fDebug)
66 {
67   //copy ctor
68 }
69
70 //____________________________________________________________________________
71 AliEMCALReconstructor::~AliEMCALReconstructor()
72 {
73   // dtor
74
75
76 //____________________________________________________________________________
77 void AliEMCALReconstructor::InitRecParam() const
78 {
79   // Check if the instance of AliEMCALRecParam exists, 
80   // if not, get it from OCDB if available, otherwise create a default one
81
82  if (!fgkRecParam  && (AliCDBManager::Instance()->IsDefaultStorageSet())) {
83     AliCDBEntry *entry = (AliCDBEntry*) 
84       AliCDBManager::Instance()->Get("EMCAL/Config/RecParam");
85     if (entry) fgkRecParam =  (AliEMCALRecParam*) entry->GetObject();
86   }
87   
88   if(!fgkRecParam){
89     AliWarning("The Reconstruction parameters for EMCAL nonitialized - Used default one");
90     fgkRecParam = new AliEMCALRecParam;
91   }
92 }
93
94 //____________________________________________________________________________
95 void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
96 {
97   // method called by AliReconstruction; 
98   // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
99   // segment maker needs access to the AliESD object to retrieve the tracks reconstructed by 
100   // the global tracking.
101   // Works on the current event.
102
103   InitRecParam();
104   AliEMCALClusterizerv1 clu;
105   clu.SetInput(digitsTree);
106   clu.SetOutput(clustersTree);
107   if(Debug())
108     clu.Digits2Clusters("deb all") ;
109   else
110     clu.Digits2Clusters("") ;
111   
112 }
113
114 //____________________________________________________________________________
115 void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
116
117 {
118   // Conversion from raw data to
119   // EMCAL digits.
120   // Works on a single-event basis
121
122   rawReader->Reset() ; 
123
124   TClonesArray *digitsArr = new TClonesArray("AliEMCALDigit",100);
125   Int_t bufsize = 32000;
126   digitsTree->Branch("EMCAL", &digitsArr, bufsize);
127
128   //Get Mapping RCU files from the AliEMCALRecParam                                                          
129   const TObjArray* maps = AliEMCALRecParam::GetMappings();
130   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
131
132   AliAltroMapping * mapping[2] ; // For the moment only 2                                                    
133   for(Int_t i = 0; i < 2; i++) {
134     mapping[i] = (AliAltroMapping*)maps->At(i);
135   }
136
137   static AliEMCALRawUtils rawUtils;
138   rawUtils.SetOption(GetOption());
139   rawUtils.Raw2Digits(rawReader,digitsArr,mapping);
140
141   digitsTree->Fill();
142   digitsArr->Delete();
143   delete digitsArr;
144
145 }
146
147 //____________________________________________________________________________
148 void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, 
149                                     AliESDEvent* esd) const
150 {
151   // Called by AliReconstruct after Reconstruct() and global tracking and vertexing 
152   // Works on the current event
153   // Creates AliESDCaloCluster from AliEMCALRecPoints 
154   // and AliESDCaloCells from AliEMCALDigits
155   // Also, fills ESD with calorimeter trigger information
156
157   //######################################################
158   //#########Calculate trigger and set trigger info###########
159   //######################################################
160  
161   AliEMCALTrigger tr ;
162   //   tr.SetPatchSize(1);//create 4x4 patches
163   tr.Trigger();
164   
165   Float_t maxAmp2x2  = tr.Get2x2MaxAmplitude();
166   Float_t maxAmpnxn  = tr.GetnxnMaxAmplitude();
167   Float_t ampOutOfPatch2x2  = tr.Get2x2AmpOutOfPatch() ;
168   Float_t ampOutOfPatchnxn  = tr.GetnxnAmpOutOfPatch() ;
169
170   AliEMCALGeometry * geom = 0;
171   AliRunLoader *runLoader = AliRunLoader::GetRunLoader();
172   if (runLoader->GetAliRun() && runLoader->GetAliRun()->GetDetector("EMCAL"))
173     geom = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
174   if (geom == 0) 
175     geom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
176
177   Int_t iSM2x2      = tr.Get2x2SuperModule();
178   Int_t iSMnxn      = tr.GetnxnSuperModule();
179   Int_t iCellPhi2x2 = tr.Get2x2CellPhi();
180   Int_t iCellPhinxn = tr.GetnxnCellPhi();
181   Int_t iCellEta2x2 = tr.Get2x2CellEta();
182   Int_t iCellEtanxn = tr.GetnxnCellEta();
183
184   AliDebug(2, Form("Trigger 2x2 max amp %f, out amp %f, SM %d, iphi %d ieta %d",  maxAmp2x2, ampOutOfPatch2x2, iSM2x2,iCellPhi2x2, iCellEta2x2));
185   AliDebug(2, Form("Trigger 4x4 max amp %f , out amp %f, SM %d, iphi %d, ieta %d",  maxAmpnxn, ampOutOfPatchnxn, iSMnxn,iCellPhinxn, iCellEtanxn));
186
187   TVector3    pos2x2(-1,-1,-1);
188   TVector3    posnxn(-1,-1,-1);
189
190   Int_t iAbsId2x2 = geom->GetAbsCellIdFromCellIndexes( iSM2x2, iCellPhi2x2, iCellEta2x2) ;
191   Int_t iAbsIdnxn = geom->GetAbsCellIdFromCellIndexes( iSMnxn, iCellPhinxn, iCellEtanxn) ;
192   geom->GetGlobal(iAbsId2x2, pos2x2);
193   geom->GetGlobal(iAbsIdnxn, posnxn);
194   
195   TArrayF triggerPosition(6);
196   triggerPosition[0] = pos2x2(0) ;   
197   triggerPosition[1] = pos2x2(1) ;   
198   triggerPosition[2] = pos2x2(2) ;  
199   triggerPosition[3] = posnxn(0) ;   
200   triggerPosition[4] = posnxn(1) ;   
201   triggerPosition[5] = posnxn(2) ;  
202
203   TArrayF triggerAmplitudes(4);
204   triggerAmplitudes[0] = maxAmp2x2 ;   
205   triggerAmplitudes[1] = ampOutOfPatch2x2 ;    
206   triggerAmplitudes[2] = maxAmpnxn ;   
207   triggerAmplitudes[3] = ampOutOfPatchnxn ;   
208
209   esd->AddEMCALTriggerPosition(triggerPosition);
210   esd->AddEMCALTriggerAmplitudes(triggerAmplitudes);
211   
212   //########################################
213   //##############Fill CaloCells###############
214   //########################################
215   
216   TClonesArray *digits = new TClonesArray("AliEMCALDigit",1000);
217   TBranch *branchdig = digitsTree->GetBranch("EMCAL");
218   if (!branchdig) { 
219     AliError("can't get the branch with the PHOS digits !");
220     return;
221   }
222   branchdig->SetAddress(&digits);
223   digitsTree->GetEvent(0);
224   Int_t nDigits = digits->GetEntries(), idignew = 0 ;
225   AliDebug(1,Form("%d digits",nDigits));
226
227   AliESDCaloCells &emcCells = *(esd->GetEMCALCells());
228   emcCells.CreateContainer(nDigits);
229   emcCells.SetType(AliESDCaloCells::kEMCALCell);
230   for (Int_t idig = 0 ; idig < nDigits ; idig++) {
231     const AliEMCALDigit * dig = (const AliEMCALDigit*)digits->At(idig);
232     if(dig->GetAmp() > 0 ){
233       emcCells.SetCell(idignew,dig->GetId(),dig->GetAmp(), dig->GetTime());   
234       idignew++;
235     }
236   }
237   emcCells.SetNumberOfCells(idignew);
238   emcCells.Sort();
239
240   //------------------------------------------------------------
241   //-----------------CLUSTERS-----------------------------
242   //------------------------------------------------------------
243   TObjArray *clusters = new TObjArray(100);
244   TBranch *branch = clustersTree->GetBranch("EMCALECARP");
245   branch->SetAddress(&clusters);
246   clustersTree->GetEvent(0);
247
248   Int_t nClusters = clusters->GetEntries(),  nClustersNew=0;
249   AliDebug(1,Form("%d clusters",nClusters));
250   esd->SetFirstEMCALCluster(esd->GetNumberOfCaloClusters()); // Put after Phos clusters 
251
252   //######################################################
253   //#######################TRACK MATCHING###############
254   //######################################################
255   //Fill list of integers, each one is index of track to which the cluster belongs.
256
257   // step 1 - initialize array of matched track indexes
258   Int_t *matchedTrack = new Int_t[nClusters];
259   for (Int_t iclus = 0; iclus < nClusters; iclus++)
260     matchedTrack[iclus] = -1;  // neg. index --> no matched track
261   
262   // step 2, change the flag for all matched clusters found in tracks
263   Int_t iemcalMatch = -1;
264   Int_t endtpc = esd->GetNumberOfTracks();
265   for (Int_t itrack = 0; itrack < endtpc; itrack++) {
266     AliESDtrack * track = esd->GetTrack(itrack) ; // retrieve track
267     iemcalMatch = track->GetEMCALcluster();
268     if(iemcalMatch >= 0) matchedTrack[iemcalMatch] = itrack;
269   } 
270
271   //########################################
272   //##############Fill CaloClusters############
273   //########################################
274
275   esd->SetNumberOfEMCALClusters(nClusters);
276   for (Int_t iClust = 0 ; iClust < nClusters ; iClust++) {
277     const AliEMCALRecPoint * clust = (const AliEMCALRecPoint*)clusters->At(iClust);
278     //if(clust->GetClusterType()== AliESDCaloCluster::kEMCALClusterv1) nRP++; else nPC++;
279     if (Debug()) clust->Print();
280     // Get information from EMCAL reconstruction points
281     Float_t xyz[3];
282     TVector3 gpos;
283     clust->GetGlobalPosition(gpos);
284     for (Int_t ixyz=0; ixyz<3; ixyz++) 
285       xyz[ixyz] = gpos[ixyz];
286     Float_t elipAxis[2];
287     clust->GetElipsAxis(elipAxis);
288     
289      //Create digits lists
290     Int_t cellMult = clust->GetMultiplicity();
291     //TArrayS digiList(digitMult);
292     Float_t *amplFloat = clust->GetEnergiesList();
293     Int_t   *digitInts = clust->GetAbsId();
294     TArrayS absIdList(cellMult);
295     //Uncomment when unfolding is done
296     //TArrayD fracList(cellMult);
297
298     Int_t newCellMult = 0; 
299     for (Int_t iCell=0; iCell<cellMult; iCell++) {
300       if (amplFloat[iCell] > 0) {
301         absIdList[newCellMult] = (UShort_t)(digitInts[iCell]);
302         //Uncomment when unfolding is done
303         //fracList[newCellMult] = amplFloat[iCell]/emcCells.GetCellAmplitude(digitInts[iCell]);
304         newCellMult++;
305       }
306     }
307     absIdList.Set(newCellMult);
308     //Uncomment when unfolding is done
309     //fracList.Set(newCellMult);
310  
311     if(newCellMult > 0) { // accept cluster if it has some digit
312       nClustersNew++;
313       //Primaries
314       Int_t  parentMult  = 0;
315       Int_t *parentList =  clust->GetParents(parentMult);
316       
317       // fills the ESDCaloCluster
318       AliESDCaloCluster * ec = new AliESDCaloCluster() ; 
319       ec->SetClusterType(AliESDCaloCluster::kEMCALClusterv1);
320       ec->SetPosition(xyz);
321       ec->SetE(clust->GetEnergy());
322       ec->SetNCells(newCellMult);
323       //Change type of list from short to ushort
324       UShort_t *newAbsIdList  = new UShort_t[newCellMult];
325       //Uncomment when unfolding is done
326       //Double_t *newFracList  = new Double_t[newCellMult];
327       for(Int_t i = 0; i < newCellMult ; i++) {
328         newAbsIdList[i]=absIdList[i];
329         //Uncomment when unfolding is done
330       //newFracList[i]=fracList[i];
331       }
332       ec->SetCellsAbsId(newAbsIdList);
333       //Uncomment when unfolding is done
334       //ec->SetCellsAmplitudeFraction(newFracList);
335       
336       ec->SetClusterDisp(clust->GetDispersion());
337       ec->SetClusterChi2(-1); //not yet implemented
338       ec->SetM02(elipAxis[0]*elipAxis[0]) ;
339       ec->SetM20(elipAxis[1]*elipAxis[1]) ;
340       ec->SetM11(-1) ;        //not yet implemented
341       
342       TArrayI arrayTrackMatched(1);// Only one track, temporal solution. 
343       arrayTrackMatched[0]= matchedTrack[iClust]; 
344       ec->AddTracksMatched(arrayTrackMatched); 
345       
346       TArrayI arrayParents(parentMult,parentList); 
347       ec->AddLabels(arrayParents);
348       
349       
350       // add the cluster to the esd object
351       esd->AddCaloCluster(ec);
352       delete ec;
353       //delete [] newAbsIdList ;
354       //delete [] newFracList ;
355     }
356   } // cycle on clusters
357   
358   delete [] matchedTrack;
359   
360   esd->SetNumberOfEMCALClusters(nClustersNew);
361   //if(nClustersNew != nClusters) 
362   //printf(" ##### nClusters %i -> new %i ##### \n", nClusters, nClustersNew );
363   
364   //Fill ESDCaloCluster with PID weights
365    AliEMCALPID *pid = new AliEMCALPID;
366    //pid->SetPrintInfo(kTRUE);
367    pid->SetReconstructor(kTRUE);
368    pid->RunPID(esd);
369    delete pid;
370   
371 }
372
373