rewrote AliMUONAlignment to - use AliMillePede2 instead of AliMillepede - use AliMill...
[u/mrichter/AliRoot.git] / MUON / AliMUONAlignmentTask.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 AliMUONAlignmentTask
20 /// AliAnalysisTask to align the MUON spectrometer.
21 /// The Task reads as input ESDs and feeds the MUONTracks to AliMUONAlignment.
22 /// The alignment itself is performed by AliMillePede2.
23 /// A OCDB entry is written with the alignment parameters.
24 ///
25 /// \author Javier Castillo, CEA/Saclay - Irfu/SPhN
26 /// \author Hugo Pereira Da Costa, CEA/Saclay - Irfu/SPhN
27 //-----------------------------------------------------------------------------
28
29 // this class header must always come first
30 #include "AliMUONAlignmentTask.h"
31
32 // local header must come before system headers
33 #include "AliAnalysisManager.h"
34 #include "AliAlignObjMatrix.h"
35 #include "AliAODInputHandler.h"
36 #include "AliAODHandler.h"
37 #include "AliAODEvent.h"
38 #include "AliAODHeader.h"
39 #include "AliCDBEntry.h"
40 #include "AliESDInputHandler.h"
41 #include "AliESDEvent.h"
42 #include "AliESDMuonTrack.h"
43 #include "AliCDBStorage.h"
44 #include "AliCDBManager.h"
45 #include "AliGRPManager.h"
46 #include "AliCDBMetaData.h"
47 #include "AliCDBId.h"
48 #include "AliGeomManager.h"
49
50 #include "AliMagF.h"
51 #include "AliMillePedeRecord.h"
52 #include "AliMpCDB.h"
53 #include "AliMUONCDB.h"
54 #include "AliMUONAlignment.h"
55 #include "AliMUONConstants.h"
56 #include "AliMUONRecoParam.h"
57 #include "AliMUONTrack.h"
58 #include "AliMUONTrackExtrap.h"
59 #include "AliMUONTrackParam.h"
60 #include "AliMUONGeometryTransformer.h"
61 #include "AliMUONESDInterface.h"
62
63 // system headers
64 #include <cassert>
65 #include <fstream>
66
67 // root headers
68 #include <TString.h>
69 #include <TError.h>
70 #include <TFile.h>
71 #include <TGraphErrors.h>
72 #include <TTree.h>
73 #include <TChain.h>
74 #include <TClonesArray.h>
75 #include <TGeoGlobalMagField.h>
76 #include <TGeoManager.h>
77 #include <Riostream.h>
78
79 ///\cond CLASSIMP
80 ClassImp(AliMUONAlignmentTask)
81 ///\endcond
82
83 //________________________________________________________________________
84 AliMUONAlignmentTask::AliMUONAlignmentTask( const char *name ):
85   AliAnalysisTaskSE(name),
86   fReadRecords( kFALSE ),
87   fWriteRecords( kFALSE ),
88   fDoAlignment( kTRUE ),
89   fMergeAlignmentCDBs( kTRUE ),
90   fForceBField( kFALSE ),
91   fBFieldOn( kFALSE ),
92   fUnbias( kFALSE ),
93   fAlign(0x0),
94   fNewAlignStorage( "local://ReAlignOCDB" ),
95   fOldGeoTransformer(0x0),
96   fNewGeoTransformer(0x0),
97   fLoadOCDBOnce(kFALSE),
98   fOCDBLoaded(kFALSE),
99   fEvent(0),
100   fTrackTot(0),
101   fTrackOk(0),
102   fRunNumberMin(0),
103   fRunNumberMax(0),
104   fRecords(0x0),
105   fRecordCount(0)
106 {
107   /// Default Constructor
108   // Define input and output slots here
109   // Input slot #0 works with a TChain
110   DefineInput(0, TChain::Class());
111
112   // initialize parameters ...
113   for(Int_t k=0;k<AliMUONAlignment::fNGlobal;k++)
114   {
115     fParameters[k]=0.;
116     fErrors[k]=0.;
117     fPulls[k]=0.;
118   }
119
120   // create alignment object
121   fAlign = new AliMUONAlignment();
122
123   // create old geometry transformer
124   fOldGeoTransformer = new AliMUONGeometryTransformer();
125
126 }
127
128 //________________________________________________________________________
129 AliMUONAlignmentTask::~AliMUONAlignmentTask()
130 {
131   /// destructor
132   delete fAlign;
133   delete fOldGeoTransformer;
134   delete fNewGeoTransformer;
135 }
136
137 //________________________________________________________________________
138 void AliMUONAlignmentTask::LocalInit()
139 {
140   /**
141   Local initialization, called once per task on the client machine
142   where the analysis train is assembled
143   **/
144
145   /* must run alignment when reading records */
146   if( fReadRecords && !fDoAlignment )
147   {
148
149     AliInfo( "Automatically setting fDoAlignment to kTRUE because fReadRecords is kTRUE" );
150     SetDoAlignment( kTRUE );
151
152   }
153
154   // print configuration
155   if( fRunNumberMin > 0 || fRunNumberMax > 0 )
156   { AliInfo( Form( "run range: %i - %i", fRunNumberMin, fRunNumberMax ) ); }
157
158   AliInfo( Form( "fReadRecords: %s", (fReadRecords ? "kTRUE":"kFALSE" ) ) );
159   AliInfo( Form( "fWriteRecords: %s", (fWriteRecords ? "kTRUE":"kFALSE" ) ) );
160   AliInfo( Form( "fDoAlignment: %s", (fDoAlignment ? "kTRUE":"kFALSE" ) ) );
161
162   if( fDoAlignment )
163   {
164     // merge alignment DB flag is irrelevant if not actually performing the alignemnt
165     AliInfo( Form( "fMergeAlignmentCDBs: %s", (fMergeAlignmentCDBs ? "kTRUE":"kFALSE" ) ) );
166   }
167
168   // storage elements
169   if( !fDefaultStorage.IsNull() ) AliInfo( Form( "fDefaultStorage: %s", fDefaultStorage.Data() ) );
170   if( !fOldAlignStorage.IsNull() ) AliInfo( Form( "fOldAlignStorage: %s", fOldAlignStorage.Data() ) );
171
172   if( fDoAlignment )
173   {
174     // new alignment storage is irrelevant if not actually performing the alignemnt
175     if( !fNewAlignStorage.IsNull() ) AliInfo( Form( "fNewAlignStorage: %s", fNewAlignStorage.Data() ) );
176     else AliFatal( "Invalid new alignment storage path" );
177   }
178
179   if( !fReadRecords )
180   {
181     // following flags are only relevant if not reading records
182     if( fForceBField ) AliInfo( Form( "fBFieldOn: %s", (fBFieldOn ? "kTRUE":"kFALSE" ) ) );
183     else AliInfo( "fBFieldOn: from GRP" );
184     AliInfo( Form( "fUnbias: %s", (fUnbias ? "kTRUE":"kFALSE" ) ) );
185   }
186
187   // consistency checks between flags
188   /* need at least one of the flags set to true */
189   if( !( fReadRecords || fWriteRecords || fDoAlignment ) )
190   { AliFatal( "Need at least one of the three flags fReadRecords, fWriteRecords and fDoAlignment set to kTRUE" ); }
191
192   /* cannot read and write records at the same time */
193   if( fReadRecords && fWriteRecords )
194   { AliFatal( "Cannot set both fReadRecords and fWriteRecords to kTRUE at the same time" ); }
195
196   /*
197   fix detectors
198   warning, counting chambers from 1.
199   this must be done before calling the Init method
200   */
201   if( fDoAlignment )
202   {
203
204     fAlign->FixChamber(1);
205     fAlign->FixChamber(10);
206
207   } else {
208
209     AliInfo( "Not fixing detector elements, since alignment is not required" );
210
211   }
212
213   // initialize
214   fAlign->Init();
215
216   // use unbiased residuals
217   fAlign->SetUnbias( fUnbias );
218
219   // Do alignment with magnetic field off
220   fAlign->SetBFieldOn( fBFieldOn );
221
222   // Set expected resolution (see AliMUONAlignment)
223   fAlign->SetSigmaXY(0.15,0.10);
224
225   // initialize global parameters to provided values
226   fAlign->InitGlobalParameters(fParameters);
227
228 }
229
230 //________________________________________________________________________
231 void AliMUONAlignmentTask::UserCreateOutputObjects()
232 {
233
234   // connect AOD output
235   if( fWriteRecords )
236   {
237     AliAODHandler* handler = dynamic_cast<AliAODHandler*>( AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler() );
238     if( handler )
239     {
240
241       // get AOD output handler and add Branch
242       fRecords = new TClonesArray( "AliMillePedeRecord", 0 );
243       fRecords->SetName( "records" );
244       handler->AddBranch("TClonesArray", &fRecords);
245
246       fRecordCount = 0;
247
248     } else AliFatal( "Error: invalid output event handler" );
249
250   }
251
252 }
253
254 //________________________________________________________________________
255 void AliMUONAlignmentTask::UserExec(Option_t *)
256 {
257
258   // do nothing if run number not in range
259   if( fRunNumberMin > 0 && fCurrentRunNumber < fRunNumberMin ) return;
260   if( fRunNumberMax > 0 && fCurrentRunNumber > fRunNumberMax ) return;
261
262   // increase event count
263   ++fEvent;
264
265   // clear array
266   if( fWriteRecords && fRecords )
267   {
268     fRecords->Clear();
269     fRecordCount = 0;
270   }
271
272   if( fReadRecords )
273   {
274
275     AliAODEvent* lAOD( dynamic_cast<AliAODEvent*>(InputEvent() ) );
276
277     // check AOD
278     if( !lAOD )
279     {
280       AliInfo("Error: AOD event not available");
281       return;
282     }
283
284     // read records
285     TClonesArray* records = static_cast<TClonesArray*>( lAOD->FindListObject( "records" ) );
286     if( !records )
287     {
288       AliInfo( "Unable to read object name \"records\"" );
289       return;
290     }
291
292     // get number of records and print
293     const Int_t lRecordsCount( records->GetEntriesFast() );
294     fTrackTot += lRecordsCount;
295
296     // loop over records
297     for( Int_t index = 0; index < lRecordsCount; ++index )
298     {
299
300       if( AliMillePedeRecord* record = dynamic_cast<AliMillePedeRecord*>( records->UncheckedAt(index) ) )
301       {
302
303         fAlign->ProcessTrack( record );
304         ++fTrackOk;
305
306         if(!(fTrackTot%100))
307         { AliInfo( Form( "Processed %i Tracks and %i were fitted.", fTrackTot, fTrackOk ) ); }
308
309       } else AliInfo( Form( "Invalid record at %i", index ) );
310
311     }
312
313   } else {
314
315     /// Main loop, called for each event
316     AliESDEvent* lESD( dynamic_cast<AliESDEvent*>(InputEvent()) );
317
318     // check ESD
319     if( !lESD )
320     {
321       AliInfo("Error: ESD event not available");
322       return;
323     }
324
325     Int_t nTracks = Int_t(lESD->GetNumberOfMuonTracks());
326     for( Int_t iTrack = 0; iTrack < nTracks; iTrack++ )
327     {
328
329       AliESDMuonTrack* esdTrack = lESD->GetMuonTrack(iTrack);
330       if (!esdTrack->ContainTrackerData()) continue;
331       if (!esdTrack->ContainTriggerData()) continue;
332
333       Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
334       if (TMath::Abs(invBenMom)<=1.04)
335       {
336
337         AliMUONTrack track;
338         AliMUONESDInterface::ESDToMUON(*esdTrack, track);
339
340         // process track and retrieve corresponding records, for storage
341         const AliMillePedeRecord* lRecords( fAlign->ProcessTrack( &track, fDoAlignment ) );
342         ++fTrackOk;
343
344         // store in array
345         if( fWriteRecords && fRecords )
346         {
347           new((*fRecords)[fRecordCount]) AliMillePedeRecord( *lRecords );
348           ++fRecordCount;
349         }
350
351       }
352
353       ++fTrackTot;
354
355       if(!(fTrackTot%100))
356       { AliInfo( Form( "Processed %i Tracks and %i were fitted.", fTrackTot, fTrackOk ) ); }
357
358       // Post final data. Write histo list to a file with option "RECREATE"
359       // PostData(1,fList);
360
361     }
362
363     // save AOD
364     if( fWriteRecords && fRecordCount > 0 )
365     {
366       AliAODHandler* handler = dynamic_cast<AliAODHandler*>( AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler() );
367       if( handler )
368       {
369         AliAODEvent* aod = handler->GetAOD();
370         AliAODHeader* header = aod->GetHeader();
371         header->SetRunNumber(lESD->GetRunNumber());
372         AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
373
374       } else AliInfo( "Error: invalid output event handler" );
375
376     }
377
378   }
379
380 }
381
382 //________________________________________________________________________
383 void AliMUONAlignmentTask::FinishTaskOutput()
384 {
385
386   /// Called once per task on the client machine at the end of the analysis.
387   AliInfo( Form( "Processed %i tracks.", fTrackTot ) );
388   AliInfo( Form( "Accepted %i tracks.", fTrackOk ) );
389
390   // stop here if no alignment is to be performed
391   if( !fDoAlignment ) return;
392
393   AliLog::SetGlobalLogLevel(AliLog::kInfo);
394
395   // Perform global fit
396   fAlign->GlobalFit(fParameters,fErrors,fPulls);
397   Int_t idOffset = 0; // 400
398   Int_t lSDetElemCh = 0;
399
400   for( Int_t iDE=0; iDE<AliMUONAlignment::fgNDetElem; iDE++ )
401   {
402
403     Int_t detElemId = idOffset+100;
404     detElemId += iDE;
405     lSDetElemCh = 0;
406     for(Int_t iCh=0; iCh<9; iCh++)
407     {
408
409       lSDetElemCh += AliMUONAlignment::fgNDetElemCh[iCh];
410       if (iDE>=lSDetElemCh)
411       {
412         detElemId += 100;
413         detElemId -= AliMUONAlignment::fgNDetElemCh[iCh];
414       }
415
416     }
417
418   }
419
420   // Post final data. Write histo list to a file with option "RECREATE"
421   // PostData(1,fList);
422
423   // store misalignments from OCDB into old transformers
424   if( fMergeAlignmentCDBs )
425   { SaveMisAlignmentData( fOldGeoTransformer ); }
426
427   // Re Align
428   fNewGeoTransformer = fAlign->ReAlign( fOldGeoTransformer, fParameters, true );
429
430   // Generate realigned data in local cdb
431   const TClonesArray* array = fNewGeoTransformer->GetMisAlignmentData();
432
433   // 100 mum residual resolution for chamber misalignments?
434   fAlign->SetAlignmentResolution( array, -1, 0.01, 0.01, 0.004, 0.003 );
435
436   // CDB manager
437   AliLog::SetGlobalDebugLevel(2);
438
439   // recover default storage full name (raw:// cannot be used to set specific storage)
440   AliCDBManager* cdbManager = AliCDBManager::Instance();
441
442   // unload old alignment path
443   if( cdbManager->GetEntryCache()->Contains("MUON/Align/Data") )
444   { cdbManager->UnloadFromCache("MUON/Align/Data"); }
445
446   // load new alignment path
447   if( !fNewAlignStorage.IsNull() ) cdbManager->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
448   else {
449
450     TString defaultStorage(cdbManager->GetDefaultStorage()->GetType());
451     if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
452     else defaultStorage += Form("://%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
453     cdbManager->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
454
455   }
456
457   // create new DB entry
458   AliCDBMetaData* cdbData = new AliCDBMetaData();
459   cdbData->SetResponsible("Dimuon Offline project");
460   cdbData->SetComment("MUON alignment objects with residual misalignment");
461   AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
462   cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
463
464 }
465
466 //________________________________________________________________________
467 void AliMUONAlignmentTask::NotifyRun()
468 {
469
470   /// run number (re)initialization
471
472   // propagate run number to fAlign
473   if( fAlign ) fAlign->SetRunNumber( fCurrentRunNumber );
474
475   // update ocdb
476   if (fOCDBLoaded && fLoadOCDBOnce)
477   {
478     AliError(Form("OCDB already loaded %d %d",fOCDBLoaded,fLoadOCDBOnce));
479     return;
480   }
481
482   AliCDBManager* cdbManager = AliCDBManager::Instance();
483
484   // Initialize default storage
485   if( !( cdbManager->IsDefaultStorageSet() || fDefaultStorage.IsNull() ) )
486   {
487
488     AliInfo( Form( "Initializing default storage: %s", fDefaultStorage.Data() ) );
489     cdbManager->SetDefaultStorage(fDefaultStorage.Data());
490
491   } else if( !fDefaultStorage.IsNull() ) {
492
493     AliInfo( "Default storage already set. Ignoring fDefaultStorage" );
494
495   } else {
496
497     AliInfo( "Default storage already set" );
498
499   }
500
501   // Initialize run number
502   if( cdbManager->GetRun() < 0 )
503   {
504     AliInfo( Form( "Setting run number: %i", fCurrentRunNumber ) );
505     cdbManager->SetRun(fCurrentRunNumber);
506   }
507
508   // following initialization is not needed when reading records
509   if( !fReadRecords )
510   {
511
512     // load magnetic field if needed
513     if( !TGeoGlobalMagField::Instance()->IsLocked() )
514     {
515
516       AliInfo( "Loading magnetic field" );
517       if( !AliMUONCDB::LoadField() )
518       {
519         AliError( "Failed to load magnetic field" );
520         return;
521       }
522
523     } else { AliInfo( "Magnetic field already locked" ); }
524
525     // checking magnitic field
526     if( !fForceBField )
527     {
528       AliInfo( "Reading magnetic field setup from GRP" );
529
530       // decide bFieldOn value base on dipole current
531       // propagete to Alignment class
532       // and printout
533       AliMagF* magF = dynamic_cast<AliMagF*>( TGeoGlobalMagField::Instance()->GetField() );
534       fBFieldOn = TMath::Abs( magF->GetFactorDip() ) > 1e-5;
535       fAlign->SetBFieldOn( fBFieldOn );
536       AliInfo( Form( "Dipole magnetic field factor: %.2f", magF->GetFactorDip() ) );
537       AliInfo( Form( "fBFieldOn = %s", (fBFieldOn ? "kTRUE":"kFALSE") ) );
538     }
539
540     AliInfo( "Loading muon mapping" );
541     if( !AliMUONCDB::LoadMapping(kFALSE) )
542     {
543       AliError( "Failed to load muon mapping" );
544       return;
545     }
546
547     AliInfo( "Assigning field to Track extrapolator" );
548     AliMUONTrackExtrap::SetField();
549
550   }
551
552   // load geometry if needed
553   if( !AliGeomManager::GetGeometry() )
554   {
555
556     // reset existing geometry/alignment if any
557     if( cdbManager->GetEntryCache()->Contains("GRP/Geometry/Data") )
558     {
559       AliInfo( "Unloading GRP/Geometry/Data" );
560       cdbManager->UnloadFromCache("GRP/Geometry/Data");
561     }
562
563     if( cdbManager->GetEntryCache()->Contains("MUON/Align/Data") )
564     {
565       AliInfo( Form( "Unloading MUON/Align/Data from %s", cdbManager->GetSpecificStorage( "MUON/Align/Data" )->GetBaseFolder().Data() ) );
566       cdbManager->UnloadFromCache("MUON/Align/Data");
567     }
568
569     // get original geometry transformer
570     AliInfo( "Loading geometry" );
571     AliGeomManager::LoadGeometry();
572     if (!AliGeomManager::GetGeometry())
573     {
574       AliError("Failed to load geometry");
575       return;
576     }
577
578     if (!fOldAlignStorage.IsNull())
579     {
580
581       AliInfo( Form( "Initializing MUON/Align/Data using: %s", fOldAlignStorage.Data() ) );
582       cdbManager->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
583
584     } else {
585
586       // recover default storage full name (raw:// cannot be used to set specific storage)
587       TString defaultStorage(cdbManager->GetDefaultStorage()->GetType());
588       if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
589       else defaultStorage += Form("://%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
590
591       AliInfo( Form( "Re-initializing MUON/Align/Data using: %s", defaultStorage.Data() ) );
592       cdbManager->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
593
594     }
595
596     AliInfo("Loading muon Alignment objects");
597     AliGeomManager::ApplyAlignObjsFromCDB("MUON");
598
599   } else { AliInfo( "Geometry already initialized - fOldAlignStorage ignored" ); }
600
601   // update geometry transformer and pass to Alignment object
602   AliInfo("Loading muon geometry in transformer");
603   fOldGeoTransformer->LoadGeometryData();
604   fAlign->SetGeometryTransformer(fOldGeoTransformer);
605
606   fOCDBLoaded = kTRUE;
607
608 }
609
610 //_____________________________________________________________________________________
611 void AliMUONAlignmentTask::SaveMisAlignmentData( AliMUONGeometryTransformer* transformer ) const
612 {
613
614   // clear transformer
615   transformer->ClearMisAlignmentData();
616
617   // load MUON/Align/Data from OCDB
618   AliCDBManager* cdbManager = AliCDBManager::Instance();
619   AliCDBEntry* cdbEntry = cdbManager->Get( "MUON/Align/Data" );
620   if (!cdbEntry)
621   {
622
623     AliError( "unable to load entry for path MUON/Align/Data" );
624     return;
625
626   }
627
628   // get TClonesArray and check
629   TClonesArray* misArray = (TClonesArray*)cdbEntry->GetObject();
630   if( !misArray )
631   {
632
633     AliError( "unable to load old misalignment array" );
634     return;
635
636   }
637
638   // loop over stored entries
639   for (Int_t index=0; index<misArray->GetEntriesFast(); ++index )
640   {
641
642     // load matrix and check
643     AliAlignObjMatrix* matrix = (AliAlignObjMatrix*) misArray->At( index );
644     if( !matrix )
645     {
646       AliError( Form( "unable to load matrix for index %i", index ) );
647       continue;
648     }
649
650     // get volume ID
651     AliGeomManager::ELayerID layerId;
652     Int_t moduleId;
653     matrix->GetVolUID( layerId, moduleId);
654
655     // make sure ELayerID is correct
656     assert( layerId == AliGeomManager::kMUON );
657
658     // printout
659     // AliInfo( Form( "Found matrix for %s %i", matrix->GetSymName(), moduleId ) );
660
661     // get matrix
662     TGeoHMatrix misMatrix;
663     matrix->GetMatrix(misMatrix);
664
665     // add to geometry transformer
666     // need the detector element
667     // "detElement->GetId()"
668     // see fOldGeoTransformer->GetMisAlignmentData( ... )
669
670     if( TString( matrix->GetSymName() ).Contains( "DE" ) ) transformer->AddMisAlignDetElement( moduleId,  misMatrix);
671     else transformer->AddMisAlignModule( moduleId,  misMatrix);
672   }
673
674   return;
675 }
676