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