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