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