fixed error of digit discarding when fill ESD cluster
[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 "AliESD.h"
32 #include "AliRunLoader.h"
33 #include "AliEMCALLoader.h"
34 #include "AliEMCALClusterizerv1.h"
35 #include "AliEMCALRecPoint.h"
36 #include "AliEMCALPID.h"
37 #include "AliRawReader.h"
38
39
40 ClassImp(AliEMCALReconstructor)
41
42 //____________________________________________________________________________
43 AliEMCALReconstructor::AliEMCALReconstructor() 
44   : fDebug(kFALSE) 
45 {
46   // ctor
47
48
49 //____________________________________________________________________________
50 AliEMCALReconstructor::AliEMCALReconstructor(const AliEMCALReconstructor & rec)
51   : AliReconstructor(rec),
52     fDebug(rec.fDebug)
53 {
54   //copy ctor
55 }
56
57 //____________________________________________________________________________
58 AliEMCALReconstructor::~AliEMCALReconstructor()
59 {
60   // dtor
61
62
63 //____________________________________________________________________________
64 void AliEMCALReconstructor::Reconstruct(AliRunLoader* runLoader) const 
65 {
66   // method called by AliReconstruction; 
67   // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
68   // segment maker needs access to the AliESD object to retrieve the tracks reconstructed by 
69   // the global tracking.
70  
71   TString headerFile(runLoader->GetFileName()) ; 
72   TString branchName(runLoader->GetEventFolder()->GetName() ) ;  
73   
74   AliEMCALClusterizerv1 clu(headerFile, branchName);
75   clu.SetEventRange(0, -1) ; // do all the events
76   if ( Debug() ) 
77     clu.ExecuteTask("deb all") ; 
78   else 
79     clu.ExecuteTask("pseudo") ;  
80  
81 }
82
83 //____________________________________________________________________________
84 void AliEMCALReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawreader) const 
85 {
86   // method called by AliReconstruction; 
87   // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
88   // segment maker needs access to the AliESD object to retrieve the tracks reconstructed by 
89   // the global tracking.
90   // Here we reconstruct from Raw Data
91   
92   rawreader->Reset() ; 
93   TString headerFile(runLoader->GetFileName()) ; 
94   TString branchName(runLoader->GetEventFolder()->GetName()) ;  
95
96   AliEMCALClusterizerv1 clu(headerFile, branchName);
97   clu.SetEventRange(0, -1) ; // do all the events
98   if ( Debug() ) 
99     clu.ExecuteTask("deb pseudo all") ; 
100   else 
101     clu.ExecuteTask("pseudo") ;  
102
103 }
104
105 //____________________________________________________________________________
106 void AliEMCALReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
107 {
108   // Called by AliReconstruct after Reconstruct() and global tracking and vertxing 
109
110   Int_t eventNumber = runLoader->GetEventNumber() ;
111
112   TString headerFile(runLoader->GetFileName()) ; 
113   TString branchName(runLoader->GetEventFolder()->GetName()) ;  
114   // Creates AliESDCaloCluster from AliEMCALRecPoints 
115   AliRunLoader *rl = AliRunLoader::GetRunLoader();
116   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
117   rl->LoadRecPoints();
118   rl->LoadKinematics(); // To get the primary label
119   rl->LoadDigits();     // To get the primary label
120   rl->LoadHits();       // To get the primary label
121   rl->GetEvent(eventNumber);
122   TObjArray *clusters = emcalLoader->RecPoints();
123   Int_t nClusters = clusters->GetEntries(), nClustersNew=0;
124   esd->SetFirstEMCALCluster(esd->GetNumberOfCaloClusters()); // Put after Phos clusters 
125   //  esd->SetNumberOfEMCALClusters(nClusters); // have to be change - Feb 25, 2007; some cluster may be discard
126
127   printf(" %i : nClusters %i \n", eventNumber, nClusters);
128   assert(0);
129   for (Int_t iClust = 0 ; iClust < nClusters ; iClust++) {
130     const AliEMCALRecPoint * clust = emcalLoader->RecPoint(iClust);
131
132     if (Debug()) clust->Print();
133     // Get information from EMCAL reconstruction points
134     Float_t xyz[3];
135     TVector3 gpos;
136     clust->GetGlobalPosition(gpos);
137     for (Int_t ixyz=0; ixyz<3; ixyz++) 
138       xyz[ixyz] = gpos[ixyz];
139
140     Int_t digitMult = clust->GetMultiplicity();
141     UShort_t *amplList = new UShort_t[digitMult];
142     UShort_t *timeList = new UShort_t[digitMult];
143     UShort_t *digiList = new UShort_t[digitMult];
144     Float_t *amplFloat = clust->GetEnergiesList();
145     Float_t *timeFloat = clust->GetTimeList();
146     Int_t   *digitInts = clust->GetAbsId();
147     Float_t elipAxis[2];
148     clust->GetElipsAxis(elipAxis);
149
150     // Convert Float_t* and Int_t* to UShort_t* to save memory
151     // Problem : we should recalculate a cluster characteristics when discard digit(s)
152     Int_t newdigitMult = 0;
153     for (Int_t iDigit=0; iDigit<digitMult; iDigit++) {
154       if(timeFloat[iDigit] < 65536/1e9*100) {
155         amplList[newdigitMult] = (UShort_t)(amplFloat[iDigit]*500);
156         if(amplList[newdigitMult] > 0) { // accept digit if poztive amplitude
157           timeList[newdigitMult] = (UShort_t)(timeFloat[iDigit]*1e9*100); // Time in units of 100 ns = 0.1 ps
158           digiList[newdigitMult] = (UShort_t)(digitInts[iDigit]);
159           newdigitMult++;
160         }
161       }
162     }
163
164     if(newdigitMult > 0) { // accept cluster if it has some digit
165       nClustersNew++;
166       if(newdigitMult != digitMult) { // some digits were deleted
167         UShort_t *amplListNew = new UShort_t[newdigitMult];
168         UShort_t *timeListNew = new UShort_t[newdigitMult];
169         UShort_t *digiListNew = new UShort_t[newdigitMult];
170         for (Int_t iDigit=0; iDigit<newdigitMult; iDigit++) {
171           amplListNew[iDigit] = amplList[iDigit];
172           timeListNew[iDigit] = timeList[iDigit];
173           digiListNew[iDigit] = digiList[iDigit];
174         }
175
176         delete [] amplList;
177         delete [] timeList;
178         delete [] digiList;
179
180         amplList = amplListNew;
181         timeList = timeListNew;
182         digiList = digiListNew;
183       }
184       // fills the ESDCaloCluster
185       AliESDCaloCluster * ec = new AliESDCaloCluster() ; 
186       ec->SetClusterType(clust->GetClusterType());
187       ec->SetGlobalPosition(xyz);
188       ec->SetClusterEnergy(clust->GetEnergy());
189
190       ec->SetNumberOfDigits(newdigitMult);
191       ec->SetDigitAmplitude(amplList); //energies
192       ec->SetDigitTime(timeList);      //times
193       ec->SetDigitIndex(digiList);     //indices
194       if(clust->GetClusterType()== AliESDCaloCluster::kClusterv1){
195         ec->SetClusterDisp(clust->GetDispersion());
196         ec->SetClusterChi2(-1); //not yet implemented
197         ec->SetM02(elipAxis[0]*elipAxis[0]) ;
198         ec->SetM20(elipAxis[1]*elipAxis[1]) ;
199         ec->SetM11(-1) ;        //not yet implemented
200         ec->SetPrimaryIndex(clust->GetPrimaryIndex());
201       } 
202     // add the cluster to the esd object
203       esd->AddCaloCluster(ec);
204       delete ec;
205     } else { // no new ESD cluster
206         delete [] amplList;
207         delete [] timeList;
208         delete [] digiList;
209     }
210   } // cycle on clusters
211   esd->SetNumberOfEMCALClusters(nClustersNew);
212   if(nClustersNew != nClusters) 
213   printf(" ##### nClusters %i -> new %i ##### \n", nClusters, nClustersNew );
214
215   //Fill ESDCaloCluster with PID weights
216   AliEMCALPID *pid = new AliEMCALPID;
217   //pid->SetPrintInfo(kTRUE);
218   pid->SetReconstructor(kTRUE);
219   pid->RunPID(esd);
220 }
221