Adding the array of the recosntruction parameters (Marian)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSReconstructor.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 "AliLog.h"
30 #include "AliAltroMapping.h"
31 #include "AliESDEvent.h"
32 #include "AliESDCaloCluster.h"
33 #include "AliESDCaloCells.h"
34 #include "AliPHOSReconstructor.h"
35 #include "AliPHOSClusterizerv1.h"
36 #include "AliPHOSTrackSegmentMakerv1.h"
37 #include "AliPHOSPIDv1.h"
38 #include "AliPHOSTracker.h"
39 #include "AliRawReader.h"
40 #include "AliCDBEntry.h"
41 #include "AliCDBManager.h"
42 #include "AliPHOSTrigger.h"
43 #include "AliPHOSGeometry.h"
44 #include "AliPHOSDigit.h"
45 #include "AliPHOSTrackSegment.h"
46 #include "AliPHOSEmcRecPoint.h"
47 #include "AliPHOSRecParticle.h"
48 #include "AliPHOSRawDecoder.h"
49 #include "AliPHOSRawDecoderv1.h"
50 #include "AliPHOSRawDecoderv2.h"
51 #include "AliPHOSRawDigiProducer.h"
52 #include "AliPHOSPulseGenerator.h"
53
54 ClassImp(AliPHOSReconstructor)
55
56 Bool_t AliPHOSReconstructor::fgDebug = kFALSE ; 
57 TClonesArray*     AliPHOSReconstructor::fgDigitsArray = 0;   // Array of PHOS digits
58 TObjArray*        AliPHOSReconstructor::fgEMCRecPoints = 0;   // Array of EMC rec.points
59
60 //____________________________________________________________________________
61 AliPHOSReconstructor::AliPHOSReconstructor() :
62   fGeom(NULL),fClusterizer(NULL),fTSM(NULL),fPID(NULL)
63 {
64   // ctor
65   fGeom        = AliPHOSGeometry::GetInstance("IHEP","");
66   fClusterizer = new AliPHOSClusterizerv1      (fGeom);
67   fTSM         = new AliPHOSTrackSegmentMakerv1(fGeom);
68   fPID         = new AliPHOSPIDv1              (fGeom);
69   fgDigitsArray = new TClonesArray("AliPHOSDigit",100);
70   fgEMCRecPoints= new TObjArray(100) ;
71 }
72
73 //____________________________________________________________________________
74   AliPHOSReconstructor::~AliPHOSReconstructor()
75 {
76   // dtor
77   delete fGeom;
78   delete fClusterizer;
79   delete fTSM;
80   delete fPID;
81   delete fgDigitsArray;
82   delete fgEMCRecPoints;
83
84
85 //____________________________________________________________________________
86 void AliPHOSReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
87 {
88   // 'single-event' local reco method called by AliReconstruction; 
89   // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
90   // segment maker needs access to the AliESDEvent object to retrieve the tracks reconstructed by 
91   // the global tracking.
92
93   fClusterizer->InitParameters();
94   fClusterizer->SetInput(digitsTree);
95   fClusterizer->SetOutput(clustersTree);
96   if ( Debug() ) 
97     fClusterizer->Digits2Clusters("deb all") ; 
98   else 
99     fClusterizer->Digits2Clusters("") ;
100 }
101
102 //____________________________________________________________________________
103 void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, 
104                                    AliESDEvent* esd) const
105 {
106   // This method produces PHOS rec-particles,
107   // then it creates AliESDtracks out of them and
108   // write tracks to the ESD
109
110
111   // do current event; the loop over events is done by AliReconstruction::Run()
112   fTSM->SetESD(esd) ; 
113   fTSM->SetInput(clustersTree);
114   if ( Debug() ) 
115     fTSM->Clusters2TrackSegments("deb all") ;
116   else 
117     fTSM->Clusters2TrackSegments("") ;
118   
119   fPID->SetEnergyCorrectionOn(GetRecoParam()->GetEMCEnergyCorrectionOn());
120   
121   fPID->SetInput(clustersTree, fTSM->GetTrackSegments()) ; 
122   fPID->SetESD(esd) ; 
123   if ( Debug() ) 
124     fPID->TrackSegments2RecParticles("deb all") ;
125   else 
126     fPID->TrackSegments2RecParticles("") ;
127
128   TClonesArray *recParticles  = fPID->GetRecParticles();
129   Int_t nOfRecParticles = recParticles->GetEntriesFast();
130   
131   esd->SetNumberOfPHOSClusters(nOfRecParticles) ; 
132   esd->SetFirstPHOSCluster(esd->GetNumberOfCaloClusters()) ;
133   
134   AliDebug(2,Form("%d rec. particles, option %s",nOfRecParticles,GetOption()));
135   
136   // Read digits array
137
138   TBranch *branch = digitsTree->GetBranch("PHOS");
139   if (!branch) { 
140     AliError("can't get the branch with the PHOS digits !");
141     return;
142   }
143   branch->SetAddress(&fgDigitsArray);
144   branch->GetEntry(0);
145
146   // Get the clusters array
147
148   TBranch *emcbranch = clustersTree->GetBranch("PHOSEmcRP");
149   if (!emcbranch) { 
150     AliError("can't get the branch with the PHOS EMC clusters !");
151     return;
152   }
153
154   emcbranch->SetAddress(&fgEMCRecPoints);
155   emcbranch->GetEntry(0);
156
157   //#########Calculate trigger and set trigger info###########
158
159   AliPHOSTrigger tr ;
160   //   tr.SetPatchSize(1);//create 4x4 patches
161   tr.SetSimulation(kFALSE);
162   tr.Trigger(fgDigitsArray);
163   
164   Float_t maxAmp2x2  = tr.Get2x2MaxAmplitude();
165   Float_t maxAmpnxn  = tr.GetnxnMaxAmplitude();
166   Float_t ampOutOfPatch2x2  = tr.Get2x2AmpOutOfPatch() ;
167   Float_t ampOutOfPatchnxn  = tr.GetnxnAmpOutOfPatch() ;
168
169   Int_t iSM2x2      = tr.Get2x2SuperModule();
170   Int_t iSMnxn      = tr.GetnxnSuperModule();
171   Int_t iCrystalPhi2x2 = tr.Get2x2CrystalPhi();
172   Int_t iCrystalPhinxn = tr.GetnxnCrystalPhi();
173   Int_t iCrystalEta2x2 = tr.Get2x2CrystalEta();
174   Int_t iCrystalEtanxn = tr.GetnxnCrystalEta();
175
176   AliDebug(2, Form("Trigger 2x2 max amp %f, out amp %f, SM %d, iphi %d ieta %d",  
177                    maxAmp2x2, ampOutOfPatch2x2, iSM2x2,iCrystalPhi2x2, iCrystalEta2x2));
178   AliDebug(2, Form("Trigger 4x4 max amp %f , out amp %f, SM %d, iphi %d, ieta %d",
179                    maxAmpnxn, ampOutOfPatchnxn, iSMnxn,iCrystalPhinxn, iCrystalEtanxn));
180
181   // Attention! PHOS modules in order to calculate AbsId need to be 1-5 not 0-4 as returns trigger.
182   Int_t iRelId2x2 []= {iSM2x2+1,0,iCrystalPhi2x2,iCrystalEta2x2};
183   Int_t iAbsId2x2 =-1;
184   Int_t iRelIdnxn []= {iSMnxn+1,0,iCrystalPhinxn,iCrystalEtanxn};
185   Int_t iAbsIdnxn =-1;
186   TVector3    pos2x2(-1,-1,-1);
187   TVector3    posnxn(-1,-1,-1);
188   fGeom->RelToAbsNumbering(iRelId2x2, iAbsId2x2);
189   fGeom->RelToAbsNumbering(iRelIdnxn, iAbsIdnxn);
190   fGeom->RelPosInAlice(iAbsId2x2, pos2x2);
191   fGeom->RelPosInAlice(iAbsIdnxn, posnxn);
192
193   TArrayF triggerPosition(6);
194   triggerPosition[0] = pos2x2(0) ;   
195   triggerPosition[1] = pos2x2(1) ;   
196   triggerPosition[2] = pos2x2(2) ;  
197   triggerPosition[3] = posnxn(0) ;   
198   triggerPosition[4] = posnxn(1) ;   
199   triggerPosition[5] = posnxn(2) ;  
200
201   TArrayF triggerAmplitudes(4);
202   triggerAmplitudes[0] = maxAmp2x2 ;   
203   triggerAmplitudes[1] = ampOutOfPatch2x2 ;    
204   triggerAmplitudes[2] = maxAmpnxn ;   
205   triggerAmplitudes[3] = ampOutOfPatchnxn ;   
206
207   //esd->SetPHOSTriggerCells(triggerPosition);
208   esd->AddPHOSTriggerPosition(triggerPosition);
209   esd->AddPHOSTriggerAmplitudes(triggerAmplitudes);
210   
211
212   //########################################
213   //############# Fill CaloCells ###########
214   //########################################
215
216   Int_t nDigits = fgDigitsArray->GetEntries();
217   Int_t idignew = 0 ;
218   AliDebug(1,Form("%d digits",nDigits));
219
220   const Int_t knEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
221   AliESDCaloCells &phsCells = *(esd->GetPHOSCells());
222   phsCells.CreateContainer(nDigits);
223   phsCells.SetType(AliESDCaloCells::kPHOSCell);
224
225   // Add to CaloCells only EMC digits with non-zero energy 
226   for (Int_t idig = 0 ; idig < nDigits ; idig++) {
227     const AliPHOSDigit * dig = (const AliPHOSDigit*)fgDigitsArray->At(idig);
228     if(dig->GetId() <= knEMC && dig->GetEnergy() > GetRecoParam()->GetEMCMinE() ){
229       //printf("i %d; id %d; amp %f; time %e\n",
230       //idignew,dig->GetId(),dig->GetEnergy(), dig->GetTime());
231       phsCells.SetCell(idignew,dig->GetId(), dig->GetEnergy(), dig->GetTime());   
232       idignew++;
233     }
234   }
235   phsCells.SetNumberOfCells(idignew);
236   phsCells.Sort();
237
238   //########################################
239   //############## Fill CaloClusters #######
240   //########################################
241
242   for (Int_t recpart = 0 ; recpart < nOfRecParticles ; recpart++) {
243     AliPHOSRecParticle  *rp    = dynamic_cast<AliPHOSRecParticle*>(recParticles->At(recpart));
244     if (Debug()) 
245       rp->Print();
246     // Get track segment and EMC rec.point associated with this rec.particle
247     AliPHOSTrackSegment *ts    = static_cast<AliPHOSTrackSegment *>(fTSM->GetTrackSegments()
248                                                                     ->At(rp->GetPHOSTSIndex()));
249
250     AliPHOSEmcRecPoint  *emcRP = static_cast<AliPHOSEmcRecPoint *>(fgEMCRecPoints->At(ts->GetEmcIndex()));
251     AliESDCaloCluster   *ec    = new AliESDCaloCluster() ; 
252     
253     Float_t xyz[3];
254     for (Int_t ixyz=0; ixyz<3; ixyz++) 
255       xyz[ixyz] = rp->GetPos()[ixyz];
256     
257     AliDebug(2,Form("Global position xyz=(%f,%f,%f)",xyz[0],xyz[1],xyz[2]));
258    
259     // Create cell lists
260
261     Int_t     cellMult   = emcRP->GetDigitsMultiplicity();
262     Int_t    *digitsList = emcRP->GetDigitsList();
263     Float_t  *rpElist    = emcRP->GetEnergiesList() ;
264     UShort_t *absIdList  = new UShort_t[cellMult];
265     Double_t *fracList   = new Double_t[cellMult];
266
267     for (Int_t iCell=0; iCell<cellMult; iCell++) {
268       AliPHOSDigit *digit = static_cast<AliPHOSDigit *>(fgDigitsArray->At(digitsList[iCell]));
269       absIdList[iCell] = (UShort_t)(digit->GetId());
270       if (digit->GetEnergy() > 0)
271         fracList[iCell] = rpElist[iCell]/digit->GetEnergy();
272       else
273         fracList[iCell] = 0;
274     }
275
276     //Primaries
277     Int_t  primMult  = 0;
278     Int_t *primList =  emcRP->GetPrimaries(primMult);
279
280     Float_t energy;
281     if (GetRecoParam()->EMCEcore2ESD())
282       energy = emcRP->GetCoreEnergy();
283     else
284       energy = rp->Energy();
285
286     // fills the ESDCaloCluster
287     ec->SetClusterType(AliESDCaloCluster::kPHOSCluster);
288     ec->SetPosition(xyz);                       //rec.point position in MARS
289     ec->SetE(energy);                           //total or core particle energy
290     ec->SetClusterDisp(emcRP->GetDispersion()); //cluster dispersion
291     ec->SetPid(rp->GetPID()) ;                  //array of particle identification
292     ec->SetM02(emcRP->GetM2x()) ;               //second moment M2x
293     ec->SetM20(emcRP->GetM2z()) ;               //second moment M2z
294     ec->SetNExMax(emcRP->GetNExMax());          //number of local maxima
295     ec->SetEmcCpvDistance(ts->GetCpvDistance("r")); //Only radius, what about separate x,z????
296     ec->SetClusterChi2(-1);                     //not yet implemented
297     ec->SetTOF(emcRP->GetTime()); //Time of flight
298
299     //Cells contributing to clusters
300     ec->SetNCells(cellMult);
301     ec->SetCellsAbsId(absIdList);
302     ec->SetCellsAmplitudeFraction(fracList);
303
304     //Distance to the nearest bad crystal
305     ec->SetDistanceToBadChannel(emcRP->GetDistanceToBadCrystal()); 
306   
307     //Array of MC indeces
308     TArrayI arrayPrim(primMult,primList);
309     ec->AddLabels(arrayPrim);
310     
311     //Matched ESD track
312     TArrayI arrayTrackMatched(1);
313     arrayTrackMatched[0]= ts->GetTrackIndex();
314     ec->AddTracksMatched(arrayTrackMatched);
315     
316     esd->AddCaloCluster(ec);
317     delete ec;   
318     delete [] fracList;
319     delete [] absIdList;
320   }
321   fgDigitsArray ->Delete();
322   fgEMCRecPoints->Delete();
323   recParticles  ->Delete();
324 }
325
326 //____________________________________________________________________________
327 AliTracker* AliPHOSReconstructor::CreateTracker() const
328 {
329   // creates the PHOS tracker
330   return new AliPHOSTracker();
331 }
332
333 //____________________________________________________________________________
334 void  AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
335 {
336   // Converts raw data to
337   // PHOS digits
338   // Works on a single-event basis
339   rawReader->Reset() ; 
340
341   AliPHOSRawDecoder * dc ;
342
343   const TObjArray* maps = AliPHOSRecoParam::GetMappings();
344   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
345
346   AliAltroMapping *mapping[4];
347   for(Int_t i = 0; i < 4; i++) {
348     mapping[i] = (AliAltroMapping*)maps->At(i);
349   }
350
351   if      (strcmp(GetRecoParam()->EMCDecoderVersion(),"v1")==0) 
352     dc=new AliPHOSRawDecoderv1(rawReader,mapping);
353   else if (strcmp(GetRecoParam()->EMCDecoderVersion(),"v2")==0) 
354     dc=new AliPHOSRawDecoderv2(rawReader,mapping);
355   else
356     dc=new AliPHOSRawDecoder(rawReader,mapping);
357
358   dc->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals());
359   
360   TClonesArray *digits = new TClonesArray("AliPHOSDigit",1);
361   digits->SetName("DIGITS");
362   Int_t bufsize = 32000;
363   digitsTree->Branch("PHOS", &digits, bufsize);
364
365   AliPHOSRawDigiProducer pr(GetRecoParam());
366   pr.MakeDigits(digits,dc);
367
368   delete dc ;
369
370   //!!!!for debug!!!
371 /*
372   Int_t modMax=-111;
373   Int_t colMax=-111;
374   Int_t rowMax=-111;
375   Float_t eMax=-333;
376   //!!!for debug!!!
377
378   Int_t relId[4];
379   for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
380     AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit);
381     if(digit->GetEnergy()>eMax) {
382       fGeom->AbsToRelNumbering(digit->GetId(),relId);
383       eMax=digit->GetEnergy();
384       modMax=relId[0];
385       rowMax=relId[2];
386       colMax=relId[3];
387     }
388   }
389
390   AliDebug(1,Form("Digit with max. energy:  modMax %d colMax %d rowMax %d  eMax %f\n\n",
391                   modMax,colMax,rowMax,eMax));
392 */
393   digitsTree->Fill();
394   digits->Delete();
395   delete digits;
396 }