Updates to AddTask macro from Chris
[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 using std::cout;
76 using std::endl;
77 ///\cond CLASSIMP   
78 ClassImp(AliMUONReAlignTask)
79 ///\endcond
80
81 //________________________________________________________________________
82 AliMUONReAlignTask::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 //________________________________________________________________________
114 AliMUONReAlignTask::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 //________________________________________________________________________________
151 AliMUONReAlignTask& AliMUONReAlignTask::operator=(const AliMUONReAlignTask& other)
152 {
153   /// Assignment
154   if(&other == this) return *this;
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 //________________________________________________________________________
177 AliMUONReAlignTask::~AliMUONReAlignTask() 
178
179   /// Destructor
180   if (fESDInterface) delete fESDInterface;
181   if (fGeoTransformer) delete fGeoTransformer;
182   if (fNewGeoTransformer) delete fNewGeoTransformer;
183 }
184
185 //________________________________________________________________________
186 void 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 //________________________________________________________________________
209 void 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 //________________________________________________________________________
229 void 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 //________________________________________________________________________
245 void 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     
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());
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         }       
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 //________________________________________________________________________
417 void 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 //-----------------------------------------------------------------------
424 void 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   }
436     
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   }
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();
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 //-----------------------------------------------------------------------
497 UInt_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 }