In AliMUONTriggerQADataMakerRec:
[u/mrichter/AliRoot.git] / MUON / AliMUONReAlignTask.cxx
CommitLineData
81f1d3ae 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
4d610fd5 18//-----------------------------------------------------------------------------
19/// \class AliMUONReAlignTask
20/// AliAnalysisTask to realign the MUON spectrometer.
21/// The Task reads as input ESDs moves the clusters of a MUONTrack acoording
22/// to the re aligned geometry taken from a misalignment file in the OCDB
23/// and refit the track. Then it produces a AliMUONClusterInfo object for each
24/// cluster. The output is a TTree of AliMUONClusterInfo
25///
26/// \author Javier Castillo, CEA/Saclay - Irfu/SPhN
27//-----------------------------------------------------------------------------
28
29#include <fstream>
30
31#include <TString.h>
32#include <TError.h>
33#include <TTree.h>
34#include <TChain.h>
35#include <TClonesArray.h>
36#include <TRandom.h>
37#include <TGeoGlobalMagField.h>
38#include <TGeoManager.h>
39#include <Riostream.h>
40
41#include "AliAnalysisTask.h"
42#include "AliAnalysisManager.h"
43#include "AliESDInputHandler.h"
44#include "AliESDEvent.h"
45#include "AliESDMuonTrack.h"
46#include "AliMagF.h"
47#include "AliCDBManager.h"
cd8521dd 48#include "AliGRPManager.h"
4d610fd5 49#include "AliGeomManager.h"
50
51#include "AliMpConstants.h"
52#include "AliMpCDB.h"
53#include "AliMpSegmentation.h"
54#include "AliMpVSegmentation.h"
55#include "AliMpPad.h"
56#include "AliMUONCalibrationData.h"
57#include "AliMUONVCalibParam.h"
58#include "AliMUONPadInfo.h"
59#include "AliMUONClusterInfo.h"
60#include "AliMUONRecoParam.h"
61#include "AliMUONESDInterface.h"
62#include "AliMUONRefitter.h"
63#include "AliMUONRecoParam.h"
64#include "AliMUONVDigit.h"
65#include "AliMUONVCluster.h"
66#include "AliMUONTrack.h"
67#include "AliMUONVTrackStore.h"
68#include "AliMUONVDigitStore.h"
69#include "AliMUONTrackParam.h"
70#include "AliMUONTrackExtrap.h"
71#include "AliMUONGeometryTransformer.h"
72
73#include "AliMUONReAlignTask.h"
74
75///\cond CLASSIMP
76ClassImp(AliMUONReAlignTask)
77///\endcond
78
79//________________________________________________________________________
80AliMUONReAlignTask::AliMUONReAlignTask(const char *name, const char *geofilename, const char *defaultocdb, const char *misalignocdb)
81 : AliAnalysisTask(name, ""),
82 fESD(0x0),
83 fClusterInfoTree(0x0),
84 fClusterInfo(0x0),
85 fESDInterface(0x0),
86 fRefitter(0x0),
87 fRecoParam(0x0),
88 fGeoFilename(geofilename),
89 fMisAlignOCDB(misalignocdb),
90 fDefaultOCDB(defaultocdb),
91 fGeoTransformer(0x0),
92 fNewGeoTransformer(0x0),
93 fGainStore(0x0),
94 fPedStore(0x0),
95 fPrintLevel(0),
96 fLastRun(-1)
97{
98 /// Default Constructor
99 // Define input and output slots here
100 // Input slot #0 works with a TChain
101 DefineInput(0, TChain::Class());
102 // Output slot #0 writes a TTree
103 DefineOutput(0, TTree::Class());
104
105 fClusterInfo = new AliMUONClusterInfo();
106 fESDInterface = new AliMUONESDInterface();
107 fGeoTransformer = new AliMUONGeometryTransformer();
108 fNewGeoTransformer = new AliMUONGeometryTransformer();
109}
110
111//________________________________________________________________________
112AliMUONReAlignTask::AliMUONReAlignTask(const AliMUONReAlignTask& obj)
113 : AliAnalysisTask(obj),
114 fESD(0x0),
115 fClusterInfoTree(0x0),
116 fClusterInfo(0x0),
117 fESDInterface(0x0),
118 fRefitter(0x0),
119 fRecoParam(0x0),
120 fGeoFilename(""),
121 fMisAlignOCDB(""),
122 fDefaultOCDB(""),
123 fGeoTransformer(0x0),
124 fNewGeoTransformer(0x0),
125 fGainStore(0x0),
126 fPedStore(0x0),
127 fPrintLevel(0),
128 fLastRun(-1)
129{
130 /// Copy constructor
131 fESD = obj.fESD;
132 fClusterInfoTree = obj.fClusterInfoTree;
133 fClusterInfo = obj.fClusterInfo;
134 fESDInterface = obj.fESDInterface;
135 fRefitter = obj.fRefitter;
136 fRecoParam = obj.fRecoParam;
137 fGeoFilename = obj.fGeoFilename;
138 fMisAlignOCDB = obj.fMisAlignOCDB;
139 fDefaultOCDB = obj.fDefaultOCDB;
140 fGeoTransformer = obj.fGeoTransformer;
141 fNewGeoTransformer = obj.fNewGeoTransformer;
142 fGainStore = obj.fGainStore;
143 fPedStore = obj.fPedStore;
144 fPrintLevel = obj.fPrintLevel;
145 fLastRun = obj.fLastRun;
146}
147
148//________________________________________________________________________________
149AliMUONReAlignTask& AliMUONReAlignTask::operator=(const AliMUONReAlignTask& other)
150{
151 /// Assignment
152 AliAnalysisTask::operator=(other);
153 fESD = other.fESD;
154 fClusterInfoTree = other.fClusterInfoTree;
155 fClusterInfo = other.fClusterInfo;
156 fESDInterface = other.fESDInterface;
157 fRefitter = other.fRefitter;
158 fRecoParam = other.fRecoParam;
159 fGeoFilename = other.fGeoFilename;
160 fMisAlignOCDB = other.fMisAlignOCDB;
161 fDefaultOCDB = other.fDefaultOCDB;
162 fGeoTransformer = other.fGeoTransformer;
163 fNewGeoTransformer = other.fNewGeoTransformer;
164 fGainStore = other.fGainStore;
165 fPedStore = other.fPedStore;
166 fPrintLevel = other.fPrintLevel;
167 fLastRun = other.fLastRun;
168
169 return *this;
170}
171
172
173//________________________________________________________________________
174AliMUONReAlignTask::~AliMUONReAlignTask()
175{
176 /// Destructor
177 if (fESDInterface) delete fESDInterface;
178 if (fGeoTransformer) delete fGeoTransformer;
179 if (fNewGeoTransformer) delete fNewGeoTransformer;
180}
181
182//________________________________________________________________________
183void AliMUONReAlignTask::LocalInit()
184{
185 /// Local initialization, called once per task on the client machine
186 /// where the analysis train is assembled
187 AliMpCDB::LoadMpSegmentation();
188
189 // prepare the refitting
190 gRandom->SetSeed(0);
191 Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
192
193 fRefitter = new AliMUONRefitter(fRecoParam);
194 fRefitter->Connect(fESDInterface);
195
196 // Original geotransformer
197 fGeoTransformer->LoadGeometryData();
198 // Apply mis alignments
199 AliGeomManager::ApplyAlignObjsFromCDB("MUON");
200 // New geotransformer
201 fNewGeoTransformer->LoadGeometryData();
202
203}
204
205//________________________________________________________________________
206void AliMUONReAlignTask::ConnectInputData(Option_t *)
207{
208 /// Connect ESD here. Called on each input data change.
209 // Connect ESD here
210 TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
211 if (!esdTree) {
212 Printf("ERROR: Could not read chain from input slot 0");
213 }
214 else {
215 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
216 if (!esdH) {
217 Printf("ERROR: Could not get ESDInputHandler");
218 }
219 else {
220 fESD = esdH->GetEvent();
221 }
222 }
223}
224
225//________________________________________________________________________
226void AliMUONReAlignTask::CreateOutputObjects()
227{
228 /// Executed once on each worker (machine actually running the analysis code)
229 //
230 // This method has to be called INSIDE the user redefined CreateOutputObjects
231 // method, before creating each object corresponding to the output containers
232 // that are to be written to a file. This need to be done in general for the big output
233 // objects that may not fit memory during processing.
234 OpenFile(0);
235
236 fClusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
237 fClusterInfoTree->Branch("clusterInfo", &fClusterInfo, 32000, 99);
238
239}
240
241//________________________________________________________________________
242void AliMUONReAlignTask::Exec(Option_t *)
243{
244 /// Main loop, called for each event
245 if (!fESD) {
246 Printf("ERROR: fESD not available");
247 return;
248 }
249
250 // Prepare Gain and Pedestal store
251 if (!fGainStore || fLastRun!=fESD->GetRunNumber()){
252 fGainStore = AliMUONCalibrationData::CreateGains(fESD->GetRunNumber());
253 fLastRun = fESD->GetRunNumber();
254 }
255 if (!fPedStore || fLastRun!=fESD->GetRunNumber()){
256 fPedStore = AliMUONCalibrationData::CreatePedestals(fESD->GetRunNumber());
257 fLastRun = fESD->GetRunNumber();
258 }
259
260 AliMUONPadInfo padInfo;
261
262 Int_t nTracks = (Int_t)fESD->GetNumberOfMuonTracks();
263 if (nTracks < 1) return;
264
265 // load the current event
266 fESDInterface->LoadEvent(*fESD);
267 AliMUONVDigitStore* digitStore = fESDInterface->GetDigits();
268
269 Double_t lX = 0.;
270 Double_t lY = 0.;
271 Double_t lZ = 0.;
272 Double_t gX = 0.;
273 Double_t gY = 0.;
274 Double_t gZ = 0.;
275 // loop over cluster to modify their position
276 AliMUONVCluster *cluster;
277 TIter next(fESDInterface->CreateClusterIterator());
278 while ((cluster = static_cast<AliMUONVCluster*>(next()))) {
279 // cout << "Original cluster" << endl;
280 // cluster->Print();
281 gX = cluster->GetX();
282 gY = cluster->GetY();
283 gZ = cluster->GetZ();
284 fGeoTransformer->Global2Local(cluster->GetDetElemId(),gX,gY,gZ,lX,lY,lZ);
285 fNewGeoTransformer->Local2Global(cluster->GetDetElemId(),lX,lY,lZ,gX,gY,gZ);
286 cluster->SetXYZ(gX,gY,gZ);
287 // cout << "Aligned cluster" << endl;
288 // cluster->Print();
289 }
290
291 // refit the tracks from digits
292 AliMUONVTrackStore* newTrackStore = fRefitter->ReconstructFromClusters();
293
294 //----------------------------------------------//
295 // ------ fill new ESD and print results ------ //
296 //----------------------------------------------//
297 // loop over the list of ESD tracks
298 TClonesArray *esdTracks = (TClonesArray*) fESD->FindListObject("MuonTracks");
299 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
300 // get the ESD track
301 AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
302
303 // skip ghost tracks (leave them unchanged in the new ESD file)
304 if (!esdTrack->ContainTrackerData()) continue;
305
306 // get the corresponding MUON track
307 AliMUONTrack* track = fESDInterface->FindTrack(esdTrack->GetUniqueID());
308
309 // Find the corresponding re-fitted MUON track
310 AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
311
312// // replace the content of the current ESD track or remove it if the refitting has failed
313// if (newTrack) {
314// Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
315// AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex, esdInterface.GetDigits());
316// } else {
317// esdTracks->Remove(esdTrack);
318// }
319
320 // print initial and re-fitted track parameters at first cluster if any
321 if (fPrintLevel>0) {
322 cout<<" ----------------track #"<<iTrack+1<<"----------------"<<endl;
323 cout<<"before refit:"<<endl;
324 AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
325 param->Print("FULL");
326 if (fPrintLevel>1) param->GetCovariances().Print();
327 if (!newTrack) continue;
328 cout<<"after refit:"<<endl;
329 param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
330 param->Print("FULL");
331 if (fPrintLevel>1) param->GetCovariances().Print();
332 cout<<" ----------------------------------------"<<endl;
333 }
334
335 // Cluster Info part
336 UInt_t muonClusterMap = BuildClusterMap(*newTrack);
337
338 // loop over clusters
339 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->First());
340 while (trackParam) {
341 fClusterInfo->Clear("C");
342
343 // fill cluster info
344 cluster = trackParam->GetClusterPtr();
345 fClusterInfo->SetRunId(fESD->GetRunNumber());
346 fClusterInfo->SetEventId(fESD->GetEventNumberInFile());
347 fClusterInfo->SetZ(cluster->GetZ());
348 fClusterInfo->SetClusterId(cluster->GetUniqueID());
349 fClusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
350 fClusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
351 fClusterInfo->SetClusterChi2(cluster->GetChi2());
352 fClusterInfo->SetClusterCharge(cluster->GetCharge());
353
354 // fill track info
355 fClusterInfo->SetTrackId(newTrack->GetUniqueID());
356 fClusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
357 fClusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetNonBendingSlope()), TMath::ATan(trackParam->GetBendingSlope()));
358 fClusterInfo->SetTrackP(trackParam->P());
359 const TMatrixD paramCov = trackParam->GetCovariances();
360 fClusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
361 fClusterInfo->SetTrackChi2(newTrack->GetNormalizedChi2());
362 fClusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());
363 fClusterInfo->SetTrackNHits(newTrack->GetNClusters());
364 fClusterInfo->SetTrackChamberHitMap(muonClusterMap);
365
366 // fill pad info if available
367 for (Int_t i=0; i<cluster->GetNDigits(); i++) {
368 AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
369 if (!digit) continue;
370
371 // pad location
372 const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
373 GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
374 AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
375
376 // calibration parameters
377 AliMUONVCalibParam* ped = fPedStore ? static_cast<AliMUONVCalibParam*>(fPedStore->FindObject(digit->DetElemId(), digit->ManuId())) : 0x0;
378 AliMUONVCalibParam* gain = fGainStore ? static_cast<AliMUONVCalibParam*>(fGainStore->FindObject(digit->DetElemId(), digit->ManuId())) : 0x0;
379 Int_t manuChannel = digit->ManuChannel();
380 Int_t planeType = 0;
381 if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
382 planeType = 1;
383 }
384
385 // fill pad info
386 padInfo.SetPadId(digit->GetUniqueID());
387 padInfo.SetPadPlaneType(planeType);
388 padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY());
389 padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY());
390 padInfo.SetPadCharge((Double_t)digit->Charge());
391 padInfo.SetPadADC(digit->ADC());
392 padInfo.SetSaturated(digit->IsSaturated());
393 padInfo.SetCalibrated(digit->IsCalibrated());
394 padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
395 padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
396 gain->ValueAsFloatFast(manuChannel,2), gain->ValueAsFloatFast(manuChannel,3));
397
398 fClusterInfo->AddPad(padInfo);
399 }
400
401 // fill cluster info tree
402 fClusterInfoTree->Fill();
403 trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
404 }
405 }
406 // free memory
407 delete newTrackStore;
408
409
410 // Post final data. Write histo list to a file with option "RECREATE"
411 PostData(0,fClusterInfoTree);
412}
413
414//________________________________________________________________________
415void AliMUONReAlignTask::Terminate(const Option_t*)
416{
417 /// Called once per task on the client machine at the end of the analysis.
418
419}
420
421//-----------------------------------------------------------------------
422void AliMUONReAlignTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
423{
424 /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
425
426 // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
427 if (!gGeoManager) {
428 AliGeomManager::LoadGeometry(geoFilename);
429 if (!gGeoManager) {
430 Error("AliMUONReAlignTask", "getting geometry from file %s failed", "generated/galice.root");
431 return;
432 }
433 }
cd8521dd 434
4d610fd5 435 // Load mapping
436 AliCDBManager* man = AliCDBManager::Instance();
437 man->SetDefaultStorage(defaultOCDB);
438 man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
439 man->Print();
440 man->SetRun(0);
441 if ( ! AliMpCDB::LoadDDLStore() ) {
442 Error("MUONRefit","Could not access mapping from OCDB !");
443 exit(-1);
444 }
cd8521dd 445
446 // set mag field
447 if (!TGeoGlobalMagField::Instance()->GetField()) {
448 printf("Loading field map...\n");
449 AliGRPManager *grpMan = new AliGRPManager();
450 grpMan->ReadGRPEntry();
451 grpMan->SetMagField();
452 delete grpMan;
453 }
454 // set the magnetic field for track extrapolations
455 AliMUONTrackExtrap::SetField();
4d610fd5 456
457 // Load initial reconstruction parameters from OCDB
458 // reconstruction parameters
459 fRecoParam = AliMUONRecoParam::GetCosmicParam();
460
461 // digit selection
462 fRecoParam->SetPadGoodnessMask(0x400BE80);
463// TString caliboption = caliboption1;
464// if ( calib == 2 ) caliboption = caliboption2;
465// fRecoParam->SetCalibrationMode("NOGAIN"caliboption.Data());
466 fRecoParam->SetCalibrationMode("NOGAIN");
467
468 // chamber resolution (incuding misalignment)
469 for (Int_t iCh=0; iCh<10; iCh++) {
470 fRecoParam->SetDefaultNonBendingReso(iCh,0.4);
471 fRecoParam->SetDefaultBendingReso(iCh,0.4);
472 }
473 fRecoParam->SetMaxNonBendingDistanceToTrack(10.);
474 fRecoParam->SetMaxBendingDistanceToTrack(10.);
475
476 // cut on (non)bending slopes
477 //fRecoParam->SetMaxNonBendingSlope(0.6);
478 //fRecoParam->SetMaxBendingSlope(0.6);
479
480 // tracking algorithm
481 // fRecoParam->MakeMoreTrackCandidates(kTRUE);
482 fRecoParam->RequestStation(0, kFALSE);
483 fRecoParam->RequestStation(2, kFALSE);
484// fRecoParam->RequestStation(3, kFALSE);
485// fRecoParam->RequestStation(4, kFALSE);
486 fRecoParam->SetSigmaCutForTracking(7.);
487 fRecoParam->ImproveTracks(kTRUE, 7.);
488 Info("MUONRefit", "\n initial recontruction parameters:");
489 fRecoParam->Print("FULL");
490
491 AliMUONESDInterface::ResetTracker(fRecoParam);
492}
493
494//-----------------------------------------------------------------------
495UInt_t AliMUONReAlignTask::BuildClusterMap(AliMUONTrack &track)
496{
497 /// Build the map of clusters in tracking chambers
498
499 UInt_t muonClusterMap = 0;
500
501 AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
502 while (trackParam) {
503
504 muonClusterMap |= BIT(trackParam->GetClusterPtr()->GetChamberId());
505
506 trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
507 }
508
509 return muonClusterMap;
510
511}