Line data Source code
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 : #include "TGeoManager.h"
26 : #include "TGeoMatrix.h"
27 :
28 : // --- Standard library ---
29 :
30 : // --- AliRoot header files ---
31 : #include "AliLog.h"
32 : #include "AliAltroMapping.h"
33 : #include "AliESDEvent.h"
34 : #include "AliESDCaloCluster.h"
35 : #include "AliESDCaloCells.h"
36 : #include "AliPHOSReconstructor.h"
37 : #include "AliPHOSClusterizerv1.h"
38 : #include "AliPHOSTrackSegmentMakerv1.h"
39 : #include "AliPHOSPIDv1.h"
40 : #include "AliPHOSTracker.h"
41 : #include "AliRawReader.h"
42 : #include "AliPHOSCalibData.h"
43 : #include "AliCDBEntry.h"
44 : #include "AliCDBManager.h"
45 : #include "AliPHOSTrigger.h"
46 : #include "AliPHOSGeometry.h"
47 : #include "AliPHOSDigit.h"
48 : #include "AliPHOSTrackSegment.h"
49 : #include "AliPHOSEmcRecPoint.h"
50 : #include "AliPHOSCpvRecPoint.h"
51 : #include "AliPHOSRecParticle.h"
52 : #include "AliPHOSRawFitterv0.h"
53 : #include "AliPHOSRawFitterv1.h"
54 : #include "AliPHOSRawFitterv2.h"
55 : #include "AliPHOSRawFitterv3.h"
56 : #include "AliPHOSRawFitterv4.h"
57 : #include "AliPHOSRawDigiProducer.h"
58 : #include "AliPHOSCpvRawDigiProducer.h"
59 : #include "AliPHOSPulseGenerator.h"
60 : #include "AliPHOSTriggerRawDigit.h"
61 : #include "AliPHOSTriggerRawDigiProducer.h"
62 : #include "AliPHOSTriggerParameters.h"
63 : #include <fstream>
64 :
65 22 : ClassImp(AliPHOSReconstructor)
66 :
67 : Bool_t AliPHOSReconstructor::fgDebug = kFALSE ;
68 : TClonesArray* AliPHOSReconstructor::fgDigitsArray = 0; // Array of PHOS digits
69 : TObjArray* AliPHOSReconstructor::fgEMCRecPoints = 0; // Array of EMC rec.points
70 : TObjArray* AliPHOSReconstructor::fgCPVRecPoints = 0; // Array of CPV rec.points
71 : AliPHOSCalibData * AliPHOSReconstructor::fgCalibData = 0 ;
72 : TClonesArray* AliPHOSReconstructor::fgTriggerDigits = 0; // Array of PHOS trigger digits
73 :
74 : //____________________________________________________________________________
75 2 : AliPHOSReconstructor::AliPHOSReconstructor() :
76 2 : fGeom(NULL),fClusterizer(NULL),fTSM(NULL),fPID(NULL),fTmpDigLG(NULL)
77 10 : {
78 : // ctor
79 4 : fGeom = AliPHOSGeometry::GetInstance("IHEP","");
80 6 : fClusterizer = new AliPHOSClusterizerv1 (fGeom);
81 6 : fTSM = new AliPHOSTrackSegmentMakerv1(fGeom);
82 6 : fPID = new AliPHOSPIDv1 (fGeom);
83 6 : fTmpDigLG = new TClonesArray("AliPHOSDigit",100);
84 6 : fgDigitsArray = new TClonesArray("AliPHOSDigit",100);
85 6 : fgEMCRecPoints = new TObjArray(100) ;
86 6 : fgCPVRecPoints = new TObjArray(100) ;
87 2 : if (!fgCalibData)
88 6 : fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
89 :
90 6 : fgTriggerDigits = new TClonesArray("AliPHOSTriggerRawDigit",100);
91 :
92 8 : AliInfo(Form("PHOS bad channel map contains %d bad channel(s).\n",
93 : fgCalibData->GetNumOfEmcBadChannels()));
94 :
95 4 : }
96 :
97 : //____________________________________________________________________________
98 : AliPHOSReconstructor::~AliPHOSReconstructor()
99 12 : {
100 : // dtor
101 : // delete fGeom; // RS: fGeom is a static object owned by AliPHOSGeometry::fGeom, don't delete
102 4 : delete fClusterizer;
103 4 : delete fTSM;
104 4 : delete fPID;
105 4 : delete fTmpDigLG;
106 4 : delete fgDigitsArray;
107 4 : delete fgEMCRecPoints;
108 4 : delete fgCPVRecPoints;
109 4 : delete fgTriggerDigits;
110 6 : }
111 :
112 : //____________________________________________________________________________
113 : void AliPHOSReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
114 : {
115 : // 'single-event' local reco method called by AliReconstruction;
116 : // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
117 : // segment maker needs access to the AliESDEvent object to retrieve the tracks reconstructed by
118 : // the global tracking.
119 :
120 16 : fClusterizer->InitParameters();
121 8 : fClusterizer->SetInput(digitsTree);
122 8 : fClusterizer->SetOutput(clustersTree);
123 16 : if ( Debug() )
124 8 : fClusterizer->Digits2Clusters("deb all") ;
125 : else
126 8 : fClusterizer->Digits2Clusters("") ;
127 8 : }
128 :
129 : //____________________________________________________________________________
130 : void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
131 : AliESDEvent* esd) const
132 : {
133 : // This method produces PHOS rec-particles,
134 : // then it creates AliESDtracks out of them and
135 : // write tracks to the ESD
136 :
137 : // do current event; the loop over events is done by AliReconstruction::Run()
138 16 : fTSM->SetESD(esd) ;
139 8 : fTSM->SetInput(clustersTree);
140 16 : if ( Debug() )
141 8 : fTSM->Clusters2TrackSegments("deb all") ;
142 : else
143 8 : fTSM->Clusters2TrackSegments("") ;
144 :
145 8 : fPID->SetInput(clustersTree, fTSM->GetTrackSegments()) ;
146 8 : fPID->SetESD(esd) ;
147 16 : if ( Debug() )
148 8 : fPID->TrackSegments2RecParticles("deb all") ;
149 : else
150 8 : fPID->TrackSegments2RecParticles("") ;
151 :
152 8 : TClonesArray *recParticles = fPID->GetRecParticles();
153 8 : Int_t nOfRecParticles = recParticles->GetEntriesFast();
154 :
155 24 : AliDebug(2,Form("%d rec. particles, option %s",nOfRecParticles,GetOption()));
156 :
157 : // Read digits array
158 :
159 8 : TBranch *branch = digitsTree->GetBranch("PHOS");
160 8 : if (!branch) {
161 0 : AliError("can't get the branch with the PHOS digits !");
162 0 : return;
163 : }
164 8 : branch->SetAddress(&fgDigitsArray);
165 8 : branch->GetEntry(0);
166 :
167 : // Get the clusters array
168 :
169 8 : TBranch *emcbranch = clustersTree->GetBranch("PHOSEmcRP");
170 8 : if (!emcbranch) {
171 0 : AliError("can't get the branch with the PHOS EMC clusters !");
172 0 : return;
173 : }
174 :
175 8 : emcbranch->SetAddress(&fgEMCRecPoints);
176 8 : emcbranch->GetEntry(0);
177 :
178 8 : TBranch *cpvbranch = clustersTree->GetBranch("PHOSCpvRP");
179 8 : if (cpvbranch) {
180 8 : cpvbranch->SetAddress(&fgCPVRecPoints);
181 8 : cpvbranch->GetEntry(0);
182 8 : }
183 :
184 : // Trigger
185 :
186 8 : TBranch *tbranch = digitsTree->GetBranch("TPHOS");
187 8 : if (tbranch) {
188 :
189 4 : tbranch->SetAddress(&fgTriggerDigits);
190 4 : tbranch->GetEntry(0);
191 :
192 8 : AliESDCaloTrigger* trgESD = esd->GetCaloTrigger("PHOS");
193 :
194 4 : if (trgESD) {
195 4 : trgESD->Allocate(fgTriggerDigits->GetEntriesFast());
196 :
197 8 : for (Int_t i = 0; i < fgTriggerDigits->GetEntriesFast(); i++) {
198 0 : AliPHOSTriggerRawDigit* tdig = (AliPHOSTriggerRawDigit*)fgTriggerDigits->At(i);
199 :
200 0 : Int_t mod,modX,modZ;
201 0 : tdig->GetModXZ(mod,modX,modZ);
202 :
203 0 : const Int_t relId[4] = {5-mod,0,modX+1,modZ+1};
204 0 : Int_t absId;
205 :
206 0 : fGeom->RelToAbsNumbering(relId,absId);
207 0 : trgESD->Add(mod,absId,tdig->GetAmp(),0.,(Int_t*)NULL,0,tdig->GetL1Threshold(),tdig->GetType());
208 0 : }
209 4 : }
210 4 : }
211 :
212 : // //#########Calculate trigger and set trigger info###########
213 :
214 : // AliPHOSTrigger tr ;
215 : // // tr.SetPatchSize(1);//create 4x4 patches
216 : // tr.SetSimulation(kFALSE);
217 : // tr.Trigger(fgDigitsArray);
218 :
219 : // Float_t maxAmp2x2 = tr.Get2x2MaxAmplitude();
220 : // Float_t maxAmpnxn = tr.GetnxnMaxAmplitude();
221 : // Float_t ampOutOfPatch2x2 = tr.Get2x2AmpOutOfPatch() ;
222 : // Float_t ampOutOfPatchnxn = tr.GetnxnAmpOutOfPatch() ;
223 :
224 : // Int_t iSM2x2 = tr.Get2x2SuperModule();
225 : // Int_t iSMnxn = tr.GetnxnSuperModule();
226 : // Int_t iCrystalPhi2x2 = tr.Get2x2CrystalPhi();
227 : // Int_t iCrystalPhinxn = tr.GetnxnCrystalPhi();
228 : // Int_t iCrystalEta2x2 = tr.Get2x2CrystalEta();
229 : // Int_t iCrystalEtanxn = tr.GetnxnCrystalEta();
230 :
231 : // AliDebug(2, Form("Trigger 2x2 max amp %f, out amp %f, SM %d, iphi %d ieta %d",
232 : // maxAmp2x2, ampOutOfPatch2x2, iSM2x2,iCrystalPhi2x2, iCrystalEta2x2));
233 : // AliDebug(2, Form("Trigger 4x4 max amp %f , out amp %f, SM %d, iphi %d, ieta %d",
234 : // maxAmpnxn, ampOutOfPatchnxn, iSMnxn,iCrystalPhinxn, iCrystalEtanxn));
235 :
236 : // // Attention! PHOS modules in order to calculate AbsId need to be 1-5 not 0-4 as returns trigger.
237 : // Int_t iRelId2x2 []= {iSM2x2+1,0,iCrystalPhi2x2,iCrystalEta2x2};
238 : // Int_t iAbsId2x2 =-1;
239 : // Int_t iRelIdnxn []= {iSMnxn+1,0,iCrystalPhinxn,iCrystalEtanxn};
240 : // Int_t iAbsIdnxn =-1;
241 : // TVector3 pos2x2(-1,-1,-1);
242 : // TVector3 posnxn(-1,-1,-1);
243 : // fGeom->RelToAbsNumbering(iRelId2x2, iAbsId2x2);
244 : // fGeom->RelToAbsNumbering(iRelIdnxn, iAbsIdnxn);
245 : // fGeom->RelPosInAlice(iAbsId2x2, pos2x2);
246 : // fGeom->RelPosInAlice(iAbsIdnxn, posnxn);
247 :
248 : // TArrayF triggerPosition(6);
249 : // triggerPosition[0] = pos2x2(0) ;
250 : // triggerPosition[1] = pos2x2(1) ;
251 : // triggerPosition[2] = pos2x2(2) ;
252 : // triggerPosition[3] = posnxn(0) ;
253 : // triggerPosition[4] = posnxn(1) ;
254 : // triggerPosition[5] = posnxn(2) ;
255 :
256 : // TArrayF triggerAmplitudes(4);
257 : // triggerAmplitudes[0] = maxAmp2x2 ;
258 : // triggerAmplitudes[1] = ampOutOfPatch2x2 ;
259 : // triggerAmplitudes[2] = maxAmpnxn ;
260 : // triggerAmplitudes[3] = ampOutOfPatchnxn ;
261 :
262 : // //esd->SetPHOSTriggerCells(triggerPosition);
263 : // esd->AddPHOSTriggerPosition(triggerPosition);
264 : // esd->AddPHOSTriggerAmplitudes(triggerAmplitudes);
265 :
266 :
267 : //########################################
268 : //############# Fill CaloCells ###########
269 : //########################################
270 8 : Int_t nDigits = fgDigitsArray->GetEntries();
271 : Int_t idignew = 0 ;
272 24 : AliDebug(1,Form("%d digits",nDigits));
273 :
274 8 : const Int_t knEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
275 8 : AliESDCaloCells &phsCells = *(esd->GetPHOSCells());
276 8 : phsCells.CreateContainer(nDigits);
277 8 : phsCells.SetType(AliESDCaloCells::kPHOSCell);
278 :
279 : // Add to CaloCells only EMC digits with non-zero energy
280 480 : for (Int_t idig = 0 ; idig < nDigits ; idig++) {
281 232 : const AliPHOSDigit * dig = (const AliPHOSDigit*)fgDigitsArray->At(idig);
282 464 : if(dig->GetId() <= knEMC &&
283 232 : Calibrate(dig->GetEnergy(),dig->GetId()) > GetRecoParam()->GetEMCMinE() ){
284 206 : Int_t primary = dig->GetPrimary(1) ;
285 412 : phsCells.SetCell(idignew,dig->GetId(), Calibrate(dig->GetEnergy(),dig->GetId()),
286 206 : CalibrateT(dig->GetTime(),dig->GetId(),dig->IsLG()),
287 206 : primary,0.,!dig->IsLG()) ;
288 :
289 206 : idignew++;
290 206 : }
291 : //Fill CPV digits
292 232 : if(dig->GetId() > knEMC &&
293 0 : Calibrate(dig->GetEnergy(),dig->GetId()) > GetRecoParam()->GetCPVMinE() ){
294 0 : Int_t primary = dig->GetPrimary(1) ;
295 : //NB! CPV cell ID can not fit Short_t <32000
296 : //So, for CPV we subtract EMC part and set negative ID
297 0 : phsCells.SetCell(idignew,-(dig->GetId()-56*64*5), Calibrate(dig->GetEnergy(),dig->GetId()),
298 : 0.,
299 : primary,0.,0) ;
300 0 : idignew++;
301 0 : }
302 : }
303 8 : phsCells.SetNumberOfCells(idignew);
304 8 : phsCells.Sort();
305 :
306 : //########################################
307 : //############## Fill CaloClusters #######
308 : //########################################
309 :
310 36 : for (Int_t recpart = 0 ; recpart < nOfRecParticles ; recpart++) {
311 10 : AliPHOSRecParticle *rp = static_cast<AliPHOSRecParticle*>(recParticles->At(recpart));
312 10 : if (Debug())
313 0 : rp->Print();
314 : // Get track segment and EMC rec.point associated with this rec.particle
315 20 : AliPHOSTrackSegment *ts = static_cast<AliPHOSTrackSegment *>(fTSM->GetTrackSegments()
316 10 : ->At(rp->GetPHOSTSIndex()));
317 :
318 10 : AliPHOSEmcRecPoint *emcRP = static_cast<AliPHOSEmcRecPoint *>(fgEMCRecPoints->At(ts->GetEmcIndex()));
319 10 : AliESDCaloCluster *ec = new AliESDCaloCluster() ;
320 :
321 10 : Float_t xyz[3];
322 80 : for (Int_t ixyz=0; ixyz<3; ixyz++)
323 60 : xyz[ixyz] = rp->GetPos()[ixyz];
324 :
325 30 : AliDebug(2,Form("Global position xyz=(%f,%f,%f)",xyz[0],xyz[1],xyz[2]));
326 :
327 : // Create cell lists
328 :
329 10 : Int_t cellMult = emcRP->GetDigitsMultiplicity();
330 10 : Int_t *digitsList = emcRP->GetDigitsList();
331 10 : Float_t *rpElist = emcRP->GetEnergiesList() ;
332 10 : UShort_t *absIdList = new UShort_t[cellMult];
333 10 : Double_t *fracList = new Double_t[cellMult];
334 :
335 344 : for (Int_t iCell=0; iCell<cellMult; iCell++) {
336 162 : AliPHOSDigit *digit = static_cast<AliPHOSDigit *>(fgDigitsArray->At(digitsList[iCell]));
337 162 : absIdList[iCell] = (UShort_t)(digit->GetId());
338 162 : if (digit->GetEnergy() > 0)
339 162 : fracList[iCell] = rpElist[iCell]/(Calibrate(digit->GetEnergy(),digit->GetId()));
340 : else
341 0 : fracList[iCell] = 0;
342 : }
343 :
344 : //Primaries
345 10 : Int_t primMult = 0;
346 10 : Int_t *primList = emcRP->GetPrimaries(primMult);
347 :
348 : Float_t energy=0.;
349 10 : if (GetRecoParam()->EMCEcore2ESD())
350 0 : energy = emcRP->GetCoreEnergy();
351 : else
352 10 : energy = rp->Energy();
353 : //Apply nonlinearity correction
354 10 : if(GetRecoParam()->GetEMCEnergyCorrectionOn())
355 10 : energy=CorrectNonlinearity(energy) ;
356 :
357 : // fills the ESDCaloCluster
358 10 : ec->SetType(AliVCluster::kPHOSNeutral);
359 10 : ec->SetPosition(xyz); //rec.point position in MARS
360 10 : ec->SetE(energy); //total or core particle energy
361 10 : ec->SetDispersion(emcRP->GetDispersion()); //cluster dispersion
362 10 : ec->SetPID(rp->GetPID()) ; //array of particle identification
363 10 : ec->SetM02(emcRP->GetM2x()) ; //second moment M2x
364 10 : ec->SetM20(emcRP->GetM2z()) ; //second moment M2z
365 10 : ec->SetNExMax(emcRP->GetNExMax()); //number of local maxima
366 10 : ec->SetEmcCpvDistance(ts->GetCpvDistance("r")); //Only radius, what about separate x,z????
367 10 : ec->SetTrackDistance(ts->GetCpvDistance("x"),ts->GetCpvDistance("z"));
368 10 : ec->SetChi2(-1); //not yet implemented
369 10 : ec->SetTOF(emcRP->GetTime()); //Time of flight - already calibrated in EMCRecPoint
370 :
371 : //Cells contributing to clusters
372 10 : ec->SetNCells(cellMult);
373 10 : ec->SetCellsAbsId(absIdList);
374 10 : ec->SetCellsAmplitudeFraction(fracList);
375 :
376 : //Distance to the nearest bad crystal
377 10 : ec->SetDistanceToBadChannel(emcRP->GetDistanceToBadCrystal());
378 :
379 : //Array of MC indeces
380 10 : TArrayI arrayPrim(primMult,primList);
381 10 : ec->AddLabels(arrayPrim);
382 :
383 : //Matched ESD track
384 10 : if(ts->GetTrackIndex()>=0){ //there is a match
385 4 : TArrayI arrayTrackMatched(1);
386 8 : arrayTrackMatched[0]= ts->GetTrackIndex();
387 4 : ec->AddTracksMatched(arrayTrackMatched);
388 :
389 4 : Int_t index = esd->AddCaloCluster(ec);
390 :
391 : //Set pointer to this cluster in ESD track
392 4 : Int_t nt=esd->GetNumberOfTracks();
393 224 : for (Int_t itr=0; itr<nt; itr++) {
394 108 : AliESDtrack *esdTrack=esd->GetTrack(itr);
395 216 : if(!esdTrack->IsPHOS())
396 100 : continue ;
397 16 : if(esdTrack->GetPHOScluster()==-recpart){ //we store negative cluster number
398 0 : esdTrack->SetPHOScluster(index) ;
399 : //no garatie that only one track matched this cluster
400 : // break ;
401 : }
402 8 : }
403 4 : }
404 : else{ //Story empty list
405 6 : TArrayI arrayTrackMatched(0);
406 6 : ec->AddTracksMatched(arrayTrackMatched);
407 6 : esd->AddCaloCluster(ec);
408 6 : }
409 :
410 20 : delete ec;
411 20 : delete [] fracList;
412 20 : delete [] absIdList;
413 10 : }
414 :
415 : //Fill CPV clusters
416 8 : Int_t nOfCPVclu=fgCPVRecPoints->GetEntriesFast() ;
417 8 : Int_t nOfTS = fTSM->GetTrackSegments()->GetEntriesFast() ;
418 16 : for (Int_t recpoint = 0 ; recpoint < nOfCPVclu ; recpoint++) {
419 0 : AliPHOSCpvRecPoint *cpvRP = static_cast<AliPHOSCpvRecPoint *>(fgCPVRecPoints->At(recpoint));
420 0 : AliESDCaloCluster *ec = new AliESDCaloCluster() ;
421 :
422 0 : Float_t xyz[3];
423 0 : TVector3 pos ;
424 0 : fGeom->GetGlobalPHOS(cpvRP, pos) ;
425 0 : for (Int_t ixyz=0; ixyz<3; ixyz++)
426 0 : xyz[ixyz] = pos[ixyz];
427 :
428 : // Create cell lists
429 0 : Int_t cellMult = cpvRP->GetDigitsMultiplicity();
430 0 : Int_t *digitsList = cpvRP->GetDigitsList();
431 0 : Float_t *rpElist = cpvRP->GetEnergiesList() ;
432 0 : UShort_t *absIdList = new UShort_t[cellMult];
433 0 : Double_t *fracList = new Double_t[cellMult];
434 :
435 0 : for (Int_t iCell=0; iCell<cellMult; iCell++) {
436 0 : AliPHOSDigit *digit = static_cast<AliPHOSDigit *>(fgDigitsArray->At(digitsList[iCell]));
437 0 : absIdList[iCell] = (UShort_t)(digit->GetId());
438 0 : if (digit->GetEnergy() > 0)
439 0 : fracList[iCell] = rpElist[iCell]/(Calibrate(digit->GetEnergy(),digit->GetId()));
440 : else
441 0 : fracList[iCell] = 0;
442 : }
443 :
444 : //Primaries
445 0 : Int_t primMult = 0;
446 0 : Int_t *primList = cpvRP->GetPrimaries(primMult);
447 :
448 0 : Float_t energy = cpvRP->GetEnergy();
449 : // fills the ESDCaloCluster
450 0 : ec->SetType(AliVCluster::kPHOSCharged);
451 0 : ec->SetPosition(xyz); //rec.point position in MARS
452 0 : ec->SetE(energy); //total or core particle energy
453 0 : ec->SetDispersion(cpvRP->GetDispersion()); //cluster dispersion
454 0 : ec->SetM02(cpvRP->GetM2x()) ; //second moment M2x
455 0 : ec->SetM20(cpvRP->GetM2z()) ; //second moment M2z
456 0 : ec->SetNExMax(cpvRP->GetNExMax()); //number of local maxima
457 0 : ec->SetChi2(-1); //not yet implemented
458 :
459 : //Cells contributing to clusters
460 0 : ec->SetNCells(cellMult);
461 0 : ec->SetCellsAbsId(absIdList);
462 0 : ec->SetCellsAmplitudeFraction(fracList);
463 :
464 : //Distance to the nearest bad crystal
465 0 : ec->SetDistanceToBadChannel(cpvRP->GetDistanceToBadCrystal());
466 :
467 : //Array of MC indeces
468 0 : TArrayI arrayPrim(primMult,primList);
469 0 : ec->AddLabels(arrayPrim);
470 :
471 : //Matched CPV-EMC clusters
472 : //CPV TrackSerments go after EMC-related TSs
473 0 : if(recpoint+nOfRecParticles<nOfTS){
474 0 : AliPHOSTrackSegment *ts = static_cast<AliPHOSTrackSegment *>(fTSM->GetTrackSegments()
475 : ->At(recpoint+nOfRecParticles));
476 0 : if(ts){
477 0 : if(ts->GetEmcIndex()>=0){
478 0 : TArrayI arrayTrackMatched(1);
479 0 : arrayTrackMatched[0]= ts->GetEmcIndex();
480 0 : ec->AddTracksMatched(arrayTrackMatched);
481 0 : ec->SetTrackDistance(ts->GetCpvDistance("x"),ts->GetCpvDistance("z"));
482 0 : }
483 : else{//No match
484 0 : TArrayI arrayTrackMatched(0);
485 0 : ec->AddTracksMatched(arrayTrackMatched);
486 0 : ec->SetTrackDistance(999.,999.);
487 0 : }
488 : }
489 0 : }
490 0 : Int_t index = esd->AddCaloCluster(ec);
491 :
492 : //Set pointer to this cluster in ESD track
493 : // Int_t nt=esd->GetNumberOfTracks();
494 : // for (Int_t itr=0; itr<nt; itr++) {
495 : // AliESDtrack *esdTrack=esd->GetTrack(itr);
496 : // if(!esdTrack->IsPHOS())
497 : // continue ;
498 : // if(esdTrack->GetPHOScluster()==-recpart){ //we store negative cluster number
499 : // esdTrack->SetPHOScluster(index) ;
500 : //no garatie that only one track matched this cluster
501 : // break ;
502 : // }
503 : // }
504 :
505 0 : delete ec;
506 0 : delete [] fracList;
507 0 : delete [] absIdList;
508 0 : }
509 :
510 :
511 8 : fgDigitsArray ->Clear("C");
512 8 : fgEMCRecPoints->Clear("C");
513 8 : recParticles ->Clear();
514 :
515 : //Store PHOS misalignment matrixes
516 8 : FillMisalMatrixes(esd) ;
517 :
518 16 : }
519 :
520 : //____________________________________________________________________________
521 : AliTracker* AliPHOSReconstructor::CreateTracker() const
522 : {
523 : // creates the PHOS tracker
524 6 : return new AliPHOSTracker();
525 0 : }
526 :
527 : //____________________________________________________________________________
528 : void AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
529 : {
530 : // Converts raw data to
531 : // PHOS digits
532 : // Works on a single-event basis
533 8 : rawReader->Reset() ;
534 :
535 : // Create a new array of PHOS and CPV digits and fill it in PHOS and CPV raw data decoders
536 8 : TClonesArray *digits = new TClonesArray("AliPHOSDigit",1);
537 4 : digits->SetName("DIGITS");
538 : Int_t bufsize = 32000;
539 4 : digitsTree->Branch("PHOS", &digits, bufsize);
540 :
541 4 : ConvertDigitsEMC(rawReader,digits);
542 4 : ConvertDigitsCPV(rawReader,digits);
543 :
544 12 : AliDebug(2,Form("Number of created digits = %d",digits->GetEntriesFast()));
545 :
546 : // Create a new array of PHOS trigger digits and fill it from raw data
547 8 : TClonesArray *tdigits = new TClonesArray("AliPHOSTriggerRawDigit",1);
548 4 : tdigits->SetName("TDIGITS");
549 4 : digitsTree->Branch("TPHOS", &tdigits, bufsize);
550 :
551 4 : rawReader->Reset();
552 4 : AliPHOSTriggerRawDigiProducer tdp(rawReader);
553 :
554 8 : AliPHOSTriggerParameters* parameters = (AliPHOSTriggerParameters*)AliPHOSRecoParam::GetTriggerParameters();
555 :
556 4 : tdp.SetTriggerParameters(parameters);
557 4 : tdp.ProcessEvent(tdigits);
558 :
559 8 : if (AliLog::GetGlobalDebugLevel() == 1) {
560 : Int_t modMax=-111;
561 : Int_t colMax=-111;
562 : Int_t rowMax=-111;
563 : Float_t eMax=-333;
564 : //!!!for debug!!!
565 :
566 0 : Int_t relId[4];
567 0 : for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
568 0 : AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit);
569 0 : if(digit->GetEnergy()>eMax) {
570 0 : fGeom->AbsToRelNumbering(digit->GetId(),relId);
571 0 : eMax=digit->GetEnergy();
572 0 : modMax=relId[0];
573 0 : rowMax=relId[2];
574 0 : colMax=relId[3];
575 0 : }
576 : }
577 :
578 0 : AliDebug(1,Form("Digit with max. energy: modMax %d colMax %d rowMax %d eMax %f\n\n",
579 : modMax,colMax,rowMax,eMax));
580 0 : }
581 :
582 4 : digitsTree->Fill();
583 :
584 4 : digits->Delete();
585 8 : delete digits;
586 :
587 4 : tdigits->Delete();
588 8 : delete tdigits;
589 4 : }
590 : //==================================================================================
591 : void AliPHOSReconstructor::ConvertDigitsEMC(AliRawReader* rawReader, TClonesArray* digits) const
592 : {
593 : // Converts CPV raw data to PHOS EMC digits
594 : // Works on a single-event basis
595 : AliPHOSRawFitterv0 * fitter ;
596 :
597 8 : const TObjArray* maps = AliPHOSRecoParam::GetMappings();
598 4 : if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
599 :
600 4 : AliAltroMapping *mapping[20];
601 168 : for(Int_t i = 0; i < 20; i++) {
602 80 : mapping[i] = (AliAltroMapping*)maps->At(i);
603 : }
604 :
605 4 : if (strcmp(GetRecoParam()->EMCFitterVersion(),"v0")==0)
606 8 : fitter=new AliPHOSRawFitterv0();
607 0 : else if (strcmp(GetRecoParam()->EMCFitterVersion(),"v1")==0)
608 0 : fitter=new AliPHOSRawFitterv1();
609 0 : else if (strcmp(GetRecoParam()->EMCFitterVersion(),"v2")==0)
610 0 : fitter=new AliPHOSRawFitterv2();
611 0 : else if (strcmp(GetRecoParam()->EMCFitterVersion(),"v3")==0)
612 0 : fitter=new AliPHOSRawFitterv3();
613 : else
614 0 : fitter=new AliPHOSRawFitterv4();
615 :
616 4 : fitter->SubtractPedestals(GetRecoParam()->EMCSubtractPedestals());
617 4 : fitter->SetAmpOffset (GetRecoParam()->GetGlobalAltroOffset());
618 4 : fitter->SetAmpThreshold (GetRecoParam()->GetGlobalAltroThreshold());
619 :
620 4 : AliPHOSRawDigiProducer rdp(rawReader,mapping);
621 :
622 8 : rdp.SetEmcMinAmp(GetRecoParam()->GetEMCRawDigitThreshold()); // in ADC
623 8 : rdp.SetCpvMinAmp(GetRecoParam()->GetCPVMinE());
624 8 : rdp.SetSampleQualityCut(GetRecoParam()->GetEMCSampleQualityCut());
625 8 : rdp.SetSubtractL1phase(GetRecoParam()->GetSubtractL1phase());
626 4 : rdp.MakeDigits(digits,fTmpDigLG,fitter);
627 :
628 8 : delete fitter ;
629 4 : }
630 : //==================================================================================
631 : void AliPHOSReconstructor::ConvertDigitsCPV(AliRawReader* rawReader, TClonesArray* digits) const
632 : {
633 : // Converts CPV raw data to PHOS CPV digits
634 : // Works on a single-event basis
635 8 : AliPHOSCpvRawDigiProducer rdp(rawReader);
636 4 : rdp.SetCalibData(fgCalibData);
637 8 : rdp.SetCpvMinAmp(GetRecoParam()->GetCPVMinE());
638 4 : rdp.SetTurbo(kTRUE);
639 4 : rdp.MakeDigits(digits);
640 :
641 4 : }
642 : //==================================================================================
643 : Float_t AliPHOSReconstructor::Calibrate(Float_t amp, Int_t absId)const{
644 : // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
645 :
646 1200 : const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
647 :
648 : //Determine rel.position of the cell absolute ID
649 600 : Int_t relId[4];
650 600 : geom->AbsToRelNumbering(absId,relId);
651 600 : Int_t module=relId[0];
652 600 : Int_t row =relId[2];
653 600 : Int_t column=relId[3];
654 600 : if(relId[1]){ //CPV
655 0 : Float_t calibration = fgCalibData->GetADCchannelCpv(module,row,column);//corrected by sevdokim
656 0 : return amp*calibration ;
657 : }
658 : else{ //EMC
659 600 : Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row);
660 600 : return amp*calibration ;
661 : }
662 600 : }
663 : //==================================================================================
664 : Float_t AliPHOSReconstructor::CalibrateT(Float_t time, Int_t absId,Bool_t isLG)const{
665 : // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
666 :
667 206 : const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
668 :
669 : //Determine rel.position of the cell absolute ID
670 206 : Int_t relId[4];
671 206 : geom->AbsToRelNumbering(absId,relId);
672 206 : Int_t module=relId[0];
673 206 : Int_t row =relId[2];
674 206 : Int_t column=relId[3];
675 206 : if(relId[1]){ //CPV
676 0 : return 0. ;
677 : }
678 : else{ //EMC
679 412 : if(isLG)
680 208 : time += fgCalibData->GetLGTimeShiftEmc(module,column,row);
681 : else
682 204 : time += fgCalibData->GetTimeShiftEmc(module,column,row);
683 206 : return time ;
684 : }
685 206 : }
686 : //==================================================================================
687 : void AliPHOSReconstructor::FillMisalMatrixes(AliESDEvent* esd)const{
688 : //Store PHOS matrixes in ESD Header
689 :
690 : //Check, if matrixes was already stored
691 46 : for(Int_t mod=0 ;mod<5; mod++){
692 16 : if(esd->GetPHOSMatrix(mod)!=0)
693 6 : return ;
694 : }
695 :
696 : //Create and store matrixes
697 2 : if(!gGeoManager){
698 0 : AliError("Can not store misal. matrixes: no gGeoManager! \n") ;
699 0 : return ;
700 : }
701 : //Note, that owner of copied marixes will be header
702 2 : char path[255] ;
703 : TGeoHMatrix * m ;
704 24 : for(Int_t mod=0; mod<5; mod++){
705 10 : snprintf(path,255,"/ALIC_1/PHOS_%d",mod+1) ; //In Geometry modules numbered 1,2,.,5 without CPV
706 10 : if (gGeoManager->CheckPath(path)){
707 8 : gGeoManager->cd(path) ;
708 8 : m = gGeoManager->GetCurrentMatrix() ;
709 16 : esd->SetPHOSMatrix(new TGeoHMatrix(*m),mod) ;
710 8 : }
711 : else{
712 2 : snprintf(path,255,"/ALIC_1/PHOC_%d",mod+1) ; //In Geometry modules numbered 1,2,3 with CPV
713 2 : if (gGeoManager->CheckPath(path)){
714 0 : gGeoManager->cd(path) ;
715 0 : m = gGeoManager->GetCurrentMatrix() ;
716 0 : esd->SetPHOSMatrix(new TGeoHMatrix(*m),mod) ;
717 0 : }
718 : else{
719 2 : snprintf(path,255,"/ALIC_1/PHOH_%d",mod+1) ; //In Geometry modules numbered 1,2,.,5
720 2 : if (gGeoManager->CheckPath(path)){
721 0 : gGeoManager->cd(path) ;
722 0 : m = gGeoManager->GetCurrentMatrix() ;
723 0 : esd->SetPHOSMatrix(new TGeoHMatrix(*m),mod) ;
724 0 : }
725 : else{
726 2 : esd->SetPHOSMatrix(NULL,mod) ;
727 : }
728 : }
729 : }
730 : }
731 :
732 10 : }
733 : //==================================================================================
734 : Float_t AliPHOSReconstructor::CorrectNonlinearity(Float_t en){
735 :
736 : //For backward compatibility, if no RecoParameters found
737 20 : if(!GetRecoParam()){
738 0 : return 0.0241+1.0504*en+0.000249*en*en ;
739 : }
740 :
741 10 : if(strcmp(GetRecoParam()->GetNonlinearityCorrectionVersion(),"NoCorrection")==0){
742 0 : return en ;
743 : }
744 10 : if(strcmp(GetRecoParam()->GetNonlinearityCorrectionVersion(),"Gustavo2005")==0){
745 10 : const Float_t *par=GetRecoParam()->GetNonlinearityParams() ;
746 10 : return par[0]+par[1]*en + par[2]*en*en ;
747 : }
748 0 : if(strcmp(GetRecoParam()->GetNonlinearityCorrectionVersion(),"Henrik2010")==0){
749 0 : const Float_t *par=GetRecoParam()->GetNonlinearityParams() ;
750 0 : return en*(par[0]+par[1]*TMath::Exp(-en*par[2]))*(1.+par[3]*TMath::Exp(-en*par[4]))*(1.+par[6]/(en*en+par[5])) ;
751 : }
752 : //For backward compatibility
753 0 : if(strcmp(GetRecoParam()->GetNonlinearityCorrectionVersion(),"")==0){
754 0 : return 0.0241+1.0504*en+0.000249*en*en ;
755 : }
756 0 : return en ;
757 10 : }
758 :
759 : void AliPHOSReconstructor::readTRUParameters(AliPHOSTriggerParameters* parameters) const
760 : {
761 : //Read trigger parameters.
762 :
763 0 : TString path(gSystem->Getenv("ALICE_ROOT"));
764 0 : path += "/PHOS/macros/Trigger/OCDB/";
765 :
766 0 : for (Int_t mod = 2; mod < 5; ++mod) { // module
767 0 : for (Int_t tru = 0; tru < 4; tru++) { // tru row
768 0 : for (Int_t branch = 0; branch < 2; branch++) { // branch
769 :
770 : // Open the Appropriate pedestal file
771 0 : TString fileName = path;
772 0 : fileName += "pedestal_m";
773 0 : fileName = fileName += mod;
774 0 : fileName+="_r";
775 0 : fileName+=tru;
776 0 : fileName+="_b";
777 0 : fileName+=branch;
778 0 : fileName+=".dat";
779 0 : std::ifstream instream;
780 0 : instream.open(fileName.Data());
781 :
782 : // Read pedestals from file
783 0 : if( ! instream.is_open() )
784 0 : Printf("E-TRUPedestals: could not open %s", fileName.Data());
785 : else
786 : {
787 0 : Int_t ped[112];
788 :
789 0 : char ch_s[36];
790 0 : char *ch_s_p = ch_s;
791 : //Int_t nlines = 0;
792 :
793 : Int_t t_ped_0 =0;
794 : Int_t t_ped_1 =0;
795 : Int_t t_ped_2 =0;
796 :
797 0 : for(Int_t n=0; n<112; n++)
798 : {
799 0 : instream.getline(ch_s_p,36);
800 0 : if (ch_s_p[23]<='9' && ch_s_p[23]>='0')
801 : {
802 0 : t_ped_0 = ch_s_p[23]-'0';
803 0 : }
804 0 : else if (ch_s_p[23]>='A' && ch_s_p[23]<='Z')
805 : {
806 0 : t_ped_0 = ch_s_p[23]-'A'+10;
807 :
808 0 : }
809 :
810 0 : if (ch_s_p[22]<='9' && ch_s_p[22]>='0')
811 : {
812 0 : t_ped_1 = ch_s_p[22]-'0';
813 0 : }
814 0 : else if (ch_s_p[22]<='Z' && ch_s_p[22]>='A')
815 : {
816 0 : t_ped_1 = ch_s_p[22]-'A'+10;
817 0 : }
818 :
819 0 : if (ch_s_p[21]<='9' && ch_s_p[21]>='0')
820 : {
821 0 : t_ped_2 = ch_s_p[21]-'0';
822 0 : }
823 0 : else if (ch_s_p[21]<='Z' && ch_s_p[21]>='A')
824 : {
825 0 : t_ped_2 = ch_s_p[21]-'A'+10;
826 0 : }
827 :
828 0 : ped[n]=t_ped_2*256+t_ped_1*16+t_ped_0;
829 :
830 :
831 : }
832 0 : for (Int_t xrow = 0; xrow < 8; xrow++){
833 0 : for (Int_t zcol = 0; zcol < 14; zcol++){
834 0 : Int_t pedestal = ped[zcol*8+xrow];
835 :
836 0 : if( pedestal < 612 && pedestal > 412 ) // resonable
837 0 : parameters->SetTRUPedestal(pedestal, mod, tru, branch, xrow, zcol);
838 : else // unresonable
839 0 : continue;
840 0 : }
841 : }
842 0 : } // Ends read of pedestals from branch from file.
843 0 : instream.close();
844 0 : }// end branch
845 : }// end tru
846 :
847 : }// end for mod
848 0 : }
849 :
850 :
851 :
852 :
|
