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