minor coverity defect: added protection for self-assignment
[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   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 //________________________________________________________________________
174 AliMUONReAlignTask::~AliMUONReAlignTask() 
175
176   /// Destructor
177   if (fESDInterface) delete fESDInterface;
178   if (fGeoTransformer) delete fGeoTransformer;
179   if (fNewGeoTransformer) delete fNewGeoTransformer;
180 }
181
182 //________________________________________________________________________
183 void 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 //________________________________________________________________________
206 void 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 //________________________________________________________________________
226 void 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 //________________________________________________________________________
242 void 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         if (ped) {
395           padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
396         } else {
397           padInfo.SetPedestal(-250.,-5.);
398         }
399         if (gain) {
400           padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
401                           gain->ValueAsIntFast(manuChannel,2), gain->ValueAsIntFast(manuChannel,3));
402         } else {
403           padInfo.SetGain(-1.,-0.1,-4095,-1);
404         }       
405         fClusterInfo->AddPad(padInfo);
406       }
407           
408       // fill cluster info tree
409       fClusterInfoTree->Fill();
410       trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
411     }
412   }
413   // free memory
414   delete newTrackStore;
415     
416   
417   // Post final data. Write histo list to a file with option "RECREATE"
418   PostData(0,fClusterInfoTree);
419 }      
420
421 //________________________________________________________________________
422 void AliMUONReAlignTask::Terminate(const Option_t*)
423 {
424   /// Called once per task on the client machine at the end of the analysis.
425
426 }
427
428 //-----------------------------------------------------------------------
429 void AliMUONReAlignTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
430 {
431   /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
432   
433   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
434   if (!gGeoManager) {
435     AliGeomManager::LoadGeometry(geoFilename);
436     if (!gGeoManager) {
437       Error("AliMUONReAlignTask", "getting geometry from file %s failed", "generated/galice.root");
438       return;
439     }
440   }
441     
442   // Load mapping
443   AliCDBManager* man = AliCDBManager::Instance();
444   man->SetDefaultStorage(defaultOCDB);
445   man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
446   man->Print();
447   man->SetRun(0);
448   if ( ! AliMpCDB::LoadDDLStore() ) {
449     Error("MUONRefit","Could not access mapping from OCDB !");
450     exit(-1);
451   }
452
453   // set mag field
454   if (!TGeoGlobalMagField::Instance()->GetField()) {
455     printf("Loading field map...\n");
456     AliGRPManager *grpMan = new AliGRPManager();
457     grpMan->ReadGRPEntry();
458     grpMan->SetMagField();
459     delete grpMan;
460   }
461   // set the magnetic field for track extrapolations
462   AliMUONTrackExtrap::SetField();
463   
464   // Load initial reconstruction parameters from OCDB
465   // reconstruction parameters
466   fRecoParam = AliMUONRecoParam::GetCosmicParam();
467   
468   // digit selection
469   fRecoParam->SetPadGoodnessMask(0x400BE80);
470 //   TString caliboption = caliboption1;
471 //   if ( calib == 2 ) caliboption = caliboption2;
472 //   fRecoParam->SetCalibrationMode("NOGAIN"caliboption.Data());
473   fRecoParam->SetCalibrationMode("NOGAIN");
474   
475   // chamber resolution (incuding misalignment)
476   for (Int_t iCh=0; iCh<10; iCh++) {
477     fRecoParam->SetDefaultNonBendingReso(iCh,0.4);
478     fRecoParam->SetDefaultBendingReso(iCh,0.4);
479   }
480   fRecoParam->SetMaxNonBendingDistanceToTrack(10.);
481   fRecoParam->SetMaxBendingDistanceToTrack(10.);
482   
483   // cut on (non)bending slopes
484   //fRecoParam->SetMaxNonBendingSlope(0.6);
485   //fRecoParam->SetMaxBendingSlope(0.6);
486   
487   // tracking algorithm
488   //  fRecoParam->MakeMoreTrackCandidates(kTRUE);
489   fRecoParam->RequestStation(0, kFALSE);
490   fRecoParam->RequestStation(2, kFALSE);
491 //   fRecoParam->RequestStation(3, kFALSE);
492 //   fRecoParam->RequestStation(4, kFALSE);
493   fRecoParam->SetSigmaCutForTracking(7.);
494   fRecoParam->ImproveTracks(kTRUE, 7.);
495   Info("MUONRefit", "\n initial recontruction parameters:");
496   fRecoParam->Print("FULL");
497   
498   AliMUONESDInterface::ResetTracker(fRecoParam);  
499 }
500
501 //-----------------------------------------------------------------------
502 UInt_t AliMUONReAlignTask::BuildClusterMap(AliMUONTrack &track)
503 {
504   /// Build the map of clusters in tracking chambers
505   
506   UInt_t muonClusterMap = 0;
507   
508   AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
509   while (trackParam) {
510     
511     muonClusterMap |= BIT(trackParam->GetClusterPtr()->GetChamberId());
512     
513     trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
514   }
515   
516   return muonClusterMap;
517   
518 }