Moved signal shape routines from AliEMCAL to separate class AliEMCALRawUtils to strea...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALReconstructor.cxx
CommitLineData
f6019cda 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 ---
f6019cda 29#include "AliEMCALReconstructor.h"
5dee926e 30
69db97f5 31#include "AliRun.h"
32#include "AliEMCAL.h"
5dee926e 33#include "AliESD.h"
34#include "AliRunLoader.h"
35#include "AliEMCALLoader.h"
f6019cda 36#include "AliEMCALClusterizerv1.h"
5dee926e 37#include "AliEMCALRecPoint.h"
dc293ae9 38#include "AliEMCALPID.h"
0964c2e9 39#include "AliEMCALTrigger.h"
1d59832c 40#include "AliRawReader.h"
41
f6019cda 42
43ClassImp(AliEMCALReconstructor)
44
45//____________________________________________________________________________
18a21c7c 46AliEMCALReconstructor::AliEMCALReconstructor()
47 : fDebug(kFALSE)
f6019cda 48{
49 // ctor
f6019cda 50}
51
0a4cb131 52//____________________________________________________________________________
18a21c7c 53AliEMCALReconstructor::AliEMCALReconstructor(const AliEMCALReconstructor & rec)
54 : AliReconstructor(rec),
55 fDebug(rec.fDebug)
0a4cb131 56{
57 //copy ctor
0a4cb131 58}
f6019cda 59
60//____________________________________________________________________________
61AliEMCALReconstructor::~AliEMCALReconstructor()
62{
63 // dtor
64}
65
66//____________________________________________________________________________
67void AliEMCALReconstructor::Reconstruct(AliRunLoader* runLoader) const
68{
69 // method called by AliReconstruction;
70 // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
71 // segment maker needs access to the AliESD object to retrieve the tracks reconstructed by
72 // the global tracking.
73
74 TString headerFile(runLoader->GetFileName()) ;
b4975670 75 TString branchName(runLoader->GetEventFolder()->GetName() ) ;
f6019cda 76
77 AliEMCALClusterizerv1 clu(headerFile, branchName);
78 clu.SetEventRange(0, -1) ; // do all the events
79 if ( Debug() )
80 clu.ExecuteTask("deb all") ;
81 else
85c60a8e 82 clu.ExecuteTask("pseudo") ;
83
f6019cda 84}
85
86//____________________________________________________________________________
5dee926e 87void AliEMCALReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawreader) const
a68156e6 88{
89 // method called by AliReconstruction;
90 // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
91 // segment maker needs access to the AliESD object to retrieve the tracks reconstructed by
92 // the global tracking.
93 // Here we reconstruct from Raw Data
94
95 rawreader->Reset() ;
96 TString headerFile(runLoader->GetFileName()) ;
97 TString branchName(runLoader->GetEventFolder()->GetName()) ;
85c60a8e 98
a68156e6 99 AliEMCALClusterizerv1 clu(headerFile, branchName);
100 clu.SetEventRange(0, -1) ; // do all the events
101 if ( Debug() )
85c60a8e 102 clu.ExecuteTask("deb pseudo all") ;
a68156e6 103 else
85c60a8e 104 clu.ExecuteTask("pseudo") ;
a68156e6 105
106}
107
108//____________________________________________________________________________
f6019cda 109void AliEMCALReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
110{
111 // Called by AliReconstruct after Reconstruct() and global tracking and vertxing
69db97f5 112 static Double_t timeScale = 1.e+11; // transition constant from sec to 0.01ns (10ps)
113
114 AliDebug(1," FillESD started ");
92da3372 115
b4975670 116 Int_t eventNumber = runLoader->GetEventNumber() ;
69db97f5 117 AliEMCALGeometry * geom = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
b4975670 118
1963b290 119 TString headerFile(runLoader->GetFileName()) ;
120 TString branchName(runLoader->GetEventFolder()->GetName()) ;
85c60a8e 121 // Creates AliESDCaloCluster from AliEMCALRecPoints
5dee926e 122 AliRunLoader *rl = AliRunLoader::GetRunLoader();
123 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
124 rl->LoadRecPoints();
92da3372 125 rl->LoadKinematics(); // To get the primary label
126 rl->LoadDigits(); // To get the primary label
127 rl->LoadHits(); // To get the primary label
5dee926e 128 rl->GetEvent(eventNumber);
129 TObjArray *clusters = emcalLoader->RecPoints();
a7a5421e 130 Int_t nClusters = clusters->GetEntries(), nClustersNew=0;
5d9e7bd7 131 // Int_t nRP=0, nPC=0; // in input
a7a5421e 132 esd->SetFirstEMCALCluster(esd->GetNumberOfCaloClusters()); // Put after Phos clusters
92da3372 133
0964c2e9 134 //#########Calculate trigger and set trigger info###########
135
136 AliEMCALTrigger tr ;
137 // tr.SetPatchSize(1);//create 4x4 patches
138 tr.Trigger();
139
140 Float_t maxAmp2x2 = tr.Get2x2MaxAmplitude();
141 Float_t maxAmpnxn = tr.GetnxnMaxAmplitude();
142 Float_t ampOutOfPatch2x2 = tr.Get2x2AmpOutOfPatch() ;
143 Float_t ampOutOfPatchnxn = tr.GetnxnAmpOutOfPatch() ;
144
0964c2e9 145 Int_t iSM2x2 = tr.Get2x2SuperModule();
146 Int_t iSMnxn = tr.GetnxnSuperModule();
147 Int_t iCellPhi2x2 = tr.Get2x2CellPhi();
148 Int_t iCellPhinxn = tr.GetnxnCellPhi();
149 Int_t iCellEta2x2 = tr.Get2x2CellEta();
150 Int_t iCellEtanxn = tr.GetnxnCellEta();
151
152 AliDebug(2, Form("Trigger 2x2 max amp %f, out amp %f, SM %d, iphi %d ieta %d", maxAmp2x2, ampOutOfPatch2x2, iSM2x2,iCellPhi2x2, iCellEta2x2));
153 AliDebug(2, Form("Trigger 4x4 max amp %f , out amp %f, SM %d, iphi %d, ieta %d", maxAmpnxn, ampOutOfPatchnxn, iSMnxn,iCellPhinxn, iCellEtanxn));
154
155 TVector3 pos2x2(-1,-1,-1);
156 TVector3 posnxn(-1,-1,-1);
157
158 Int_t iAbsId2x2 = geom->GetAbsCellIdFromCellIndexes( iSM2x2, iCellPhi2x2, iCellEta2x2) ;
159 Int_t iAbsIdnxn = geom->GetAbsCellIdFromCellIndexes( iSMnxn, iCellPhinxn, iCellEtanxn) ;
160 geom->GetGlobal(iAbsId2x2, pos2x2);
161 geom->GetGlobal(iAbsIdnxn, posnxn);
162
163 TArrayF triggerPosition(6);
164 triggerPosition[0] = pos2x2(0) ;
165 triggerPosition[1] = pos2x2(1) ;
166 triggerPosition[2] = pos2x2(2) ;
167 triggerPosition[3] = posnxn(0) ;
168 triggerPosition[4] = posnxn(1) ;
169 triggerPosition[5] = posnxn(2) ;
170
171 TArrayF triggerAmplitudes(4);
172 triggerAmplitudes[0] = maxAmp2x2 ;
173 triggerAmplitudes[1] = ampOutOfPatch2x2 ;
174 triggerAmplitudes[2] = maxAmpnxn ;
175 triggerAmplitudes[3] = ampOutOfPatchnxn ;
176
177 esd->AddEMCALTriggerPosition(triggerPosition);
178 esd->AddEMCALTriggerAmplitudes(triggerAmplitudes);
179
180 //######################################
181
182 //Fill CaloClusters
183
5dee926e 184 for (Int_t iClust = 0 ; iClust < nClusters ; iClust++) {
185 const AliEMCALRecPoint * clust = emcalLoader->RecPoint(iClust);
5d9e7bd7 186 //if(clust->GetClusterType()== AliESDCaloCluster::kClusterv1) nRP++; else nPC++;
85c60a8e 187 if (Debug()) clust->Print();
a7a5421e 188 // Get information from EMCAL reconstruction points
85c60a8e 189 Float_t xyz[3];
5dee926e 190 TVector3 gpos;
191 clust->GetGlobalPosition(gpos);
f6019cda 192 for (Int_t ixyz=0; ixyz<3; ixyz++)
5dee926e 193 xyz[ixyz] = gpos[ixyz];
92da3372 194
85c60a8e 195 Int_t digitMult = clust->GetMultiplicity();
196 UShort_t *amplList = new UShort_t[digitMult];
197 UShort_t *timeList = new UShort_t[digitMult];
198 UShort_t *digiList = new UShort_t[digitMult];
199 Float_t *amplFloat = clust->GetEnergiesList();
200 Float_t *timeFloat = clust->GetTimeList();
201 Int_t *digitInts = clust->GetAbsId();
8a11aef1 202 Float_t elipAxis[2];
92da3372 203 clust->GetElipsAxis(elipAxis);
85c60a8e 204
a7a5421e 205 // Convert Float_t* and Int_t* to UShort_t* to save memory
206 // Problem : we should recalculate a cluster characteristics when discard digit(s)
f1487f22 207 Int_t newdigitMult = 0;
92da3372 208 for (Int_t iDigit=0; iDigit<digitMult; iDigit++) {
f1487f22 209 if(timeFloat[iDigit] < 65536./timeScale) {
a7a5421e 210 amplList[newdigitMult] = (UShort_t)(amplFloat[iDigit]*500);
211 if(amplList[newdigitMult] > 0) { // accept digit if poztive amplitude
f1487f22 212 timeList[newdigitMult] = (UShort_t)(timeFloat[iDigit]*timeScale); // Time in units of 0.01 ns = 10 ps
a7a5421e 213 digiList[newdigitMult] = (UShort_t)(digitInts[iDigit]);
214 newdigitMult++;
215 }
92da3372 216 }
92da3372 217 }
a7a5421e 218
219 if(newdigitMult > 0) { // accept cluster if it has some digit
220 nClustersNew++;
221 if(newdigitMult != digitMult) { // some digits were deleted
222 UShort_t *amplListNew = new UShort_t[newdigitMult];
223 UShort_t *timeListNew = new UShort_t[newdigitMult];
224 UShort_t *digiListNew = new UShort_t[newdigitMult];
225 for (Int_t iDigit=0; iDigit<newdigitMult; iDigit++) {
226 amplListNew[iDigit] = amplList[iDigit];
227 timeListNew[iDigit] = timeList[iDigit];
228 digiListNew[iDigit] = digiList[iDigit];
229 }
230
231 delete [] amplList;
232 delete [] timeList;
233 delete [] digiList;
234
235 amplList = amplListNew;
236 timeList = timeListNew;
237 digiList = digiListNew;
238 }
65721814 239
240 //Primaries
241 Int_t primMult = 0;
242 Int_t *primInts = clust->GetPrimaries(primMult);
243 UShort_t *primList = new UShort_t[primMult];
244 for (Int_t ipr=0; ipr<primMult; ipr++)
245 primList[ipr] = (UShort_t)(primInts[ipr]);
246
247
a7a5421e 248 // fills the ESDCaloCluster
249 AliESDCaloCluster * ec = new AliESDCaloCluster() ;
69db97f5 250 ec->SetEMCAL(kTRUE);
a7a5421e 251 ec->SetClusterType(clust->GetClusterType());
69db97f5 252
a7a5421e 253 ec->SetNumberOfDigits(newdigitMult);
254 ec->SetDigitAmplitude(amplList); //energies
255 ec->SetDigitTime(timeList); //times
256 ec->SetDigitIndex(digiList); //indices
257 if(clust->GetClusterType()== AliESDCaloCluster::kClusterv1){
69db97f5 258 ec->SetClusterEnergy(clust->GetEnergy());
259 ec->SetGlobalPosition(xyz);
260
65721814 261 ec->SetPrimaryIndex(clust->GetPrimaryIndex());
262 ec->SetNumberOfPrimaries(primMult); //primary multiplicity
263 ec->SetListOfPrimaries(primList); //primary List for a cluster
69db97f5 264
a7a5421e 265 ec->SetClusterDisp(clust->GetDispersion());
266 ec->SetClusterChi2(-1); //not yet implemented
267 ec->SetM02(elipAxis[0]*elipAxis[0]) ;
268 ec->SetM20(elipAxis[1]*elipAxis[1]) ;
269 ec->SetM11(-1) ; //not yet implemented
a7a5421e 270 }
85c60a8e 271 // add the cluster to the esd object
a7a5421e 272 esd->AddCaloCluster(ec);
273 delete ec;
274 } else { // no new ESD cluster
275 delete [] amplList;
276 delete [] timeList;
277 delete [] digiList;
278 }
279 } // cycle on clusters
d948db41 280 // printf(" %i : nClusters %i : RP %i PC %i \n", eventNumber, nClusters, nRP, nPC);
a7a5421e 281 esd->SetNumberOfEMCALClusters(nClustersNew);
5d9e7bd7 282 //if(nClustersNew != nClusters)
283 //printf(" ##### nClusters %i -> new %i ##### \n", nClusters, nClustersNew );
dc293ae9 284
dc293ae9 285 //Fill ESDCaloCluster with PID weights
dc293ae9 286 AliEMCALPID *pid = new AliEMCALPID;
287 //pid->SetPrintInfo(kTRUE);
288 pid->SetReconstructor(kTRUE);
289 pid->RunPID(esd);
69db97f5 290
291 AliDebug(1," FillESD ended ");
f6019cda 292}
dc293ae9 293