fill information about type of cluster on ESD - SetEMCAL(kTRUE); correct the getting...
[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 "AliRun.h"
32 #include "AliEMCAL.h"
33 #include "AliESD.h"
34 #include "AliRunLoader.h"
35 #include "AliEMCALLoader.h"
36 #include "AliEMCALClusterizerv1.h"
37 #include "AliEMCALRecPoint.h"
38 #include "AliEMCALPID.h"
39 #include "AliEMCALTrigger.h"
40 #include "AliRawReader.h"
41
42
43 ClassImp(AliEMCALReconstructor)
44
45 //____________________________________________________________________________
46 AliEMCALReconstructor::AliEMCALReconstructor() 
47   : fDebug(kFALSE) 
48 {
49   // ctor
50
51
52 //____________________________________________________________________________
53 AliEMCALReconstructor::AliEMCALReconstructor(const AliEMCALReconstructor & rec)
54   : AliReconstructor(rec),
55     fDebug(rec.fDebug)
56 {
57   //copy ctor
58 }
59
60 //____________________________________________________________________________
61 AliEMCALReconstructor::~AliEMCALReconstructor()
62 {
63   // dtor
64
65
66 //____________________________________________________________________________
67 void 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()) ; 
75   TString branchName(runLoader->GetEventFolder()->GetName() ) ;  
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 
82     clu.ExecuteTask("pseudo") ;  
83  
84 }
85
86 //____________________________________________________________________________
87 void AliEMCALReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawreader) const 
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()) ;  
98
99   AliEMCALClusterizerv1 clu(headerFile, branchName);
100   clu.SetEventRange(0, -1) ; // do all the events
101   if ( Debug() ) 
102     clu.ExecuteTask("deb pseudo all") ; 
103   else 
104     clu.ExecuteTask("pseudo") ;  
105
106 }
107
108 //____________________________________________________________________________
109 void AliEMCALReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
110 {
111   // Called by AliReconstruct after Reconstruct() and global tracking and vertxing 
112   static  Double_t timeScale = 1.e+11; // transition constant from sec to 0.01ns (10ps)
113
114   AliDebug(1," FillESD started ");
115
116   Int_t eventNumber = runLoader->GetEventNumber() ;
117   AliEMCALGeometry * geom = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
118
119   TString headerFile(runLoader->GetFileName()) ; 
120   TString branchName(runLoader->GetEventFolder()->GetName()) ;  
121   // Creates AliESDCaloCluster from AliEMCALRecPoints 
122   AliRunLoader *rl = AliRunLoader::GetRunLoader();
123   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
124   rl->LoadRecPoints();
125   rl->LoadKinematics(); // To get the primary label
126   rl->LoadDigits();     // To get the primary label
127   rl->LoadHits();       // To get the primary label
128   rl->GetEvent(eventNumber);
129   TObjArray *clusters = emcalLoader->RecPoints();
130   Int_t nClusters = clusters->GetEntries(), nClustersNew=0;
131   //  Int_t nRP=0, nPC=0; // in input
132   esd->SetFirstEMCALCluster(esd->GetNumberOfCaloClusters()); // Put after Phos clusters 
133
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
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
184   for (Int_t iClust = 0 ; iClust < nClusters ; iClust++) {
185     const AliEMCALRecPoint * clust = emcalLoader->RecPoint(iClust);
186     //if(clust->GetClusterType()== AliESDCaloCluster::kClusterv1) nRP++; else nPC++;
187     if (Debug()) clust->Print();
188     // Get information from EMCAL reconstruction points
189     Float_t xyz[3];
190     TVector3 gpos;
191     clust->GetGlobalPosition(gpos);
192     for (Int_t ixyz=0; ixyz<3; ixyz++) 
193       xyz[ixyz] = gpos[ixyz];
194
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();
202     Float_t elipAxis[2];
203     clust->GetElipsAxis(elipAxis);
204
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)
207     Int_t newdigitMult = 0; 
208     for (Int_t iDigit=0; iDigit<digitMult; iDigit++) {
209       if(timeFloat[iDigit] < 65536./timeScale) {
210         amplList[newdigitMult] = (UShort_t)(amplFloat[iDigit]*500);
211         if(amplList[newdigitMult] > 0) { // accept digit if poztive amplitude
212           timeList[newdigitMult] = (UShort_t)(timeFloat[iDigit]*timeScale); // Time in units of 0.01 ns = 10 ps
213           digiList[newdigitMult] = (UShort_t)(digitInts[iDigit]);
214           newdigitMult++;
215         }
216       }
217     }
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       }
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
248       // fills the ESDCaloCluster
249       AliESDCaloCluster * ec = new AliESDCaloCluster() ; 
250       ec->SetEMCAL(kTRUE);
251       ec->SetClusterType(clust->GetClusterType());
252
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){
258         ec->SetClusterEnergy(clust->GetEnergy());
259         ec->SetGlobalPosition(xyz);
260
261         ec->SetPrimaryIndex(clust->GetPrimaryIndex());
262         ec->SetNumberOfPrimaries(primMult);           //primary multiplicity
263         ec->SetListOfPrimaries(primList);                  //primary List for a cluster  
264
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
270       } 
271     // add the cluster to the esd object
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
280   //  printf(" %i : nClusters %i : RP %i PC %i \n", eventNumber, nClusters, nRP, nPC);
281   esd->SetNumberOfEMCALClusters(nClustersNew);
282   //if(nClustersNew != nClusters) 
283   //printf(" ##### nClusters %i -> new %i ##### \n", nClusters, nClustersNew );
284
285   //Fill ESDCaloCluster with PID weights
286   AliEMCALPID *pid = new AliEMCALPID;
287   //pid->SetPrintInfo(kTRUE);
288   pid->SetReconstructor(kTRUE);
289   pid->RunPID(esd);
290
291   AliDebug(1," FillESD ended ");
292 }
293