remove using unbiased residuals: it does not work fixed compilation warnings
[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 = aod->GetHeader();
368         header->SetRunNumber(lESD->GetRunNumber());
369         AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
370
371       } else AliInfo( "Error: invalid output event handler" );
372
373     }
374
375   }
376
377 }
378
379 //________________________________________________________________________
380 void AliMUONAlignmentTask::FinishTaskOutput()
381 {
382
383   /// Called once per task on the client machine at the end of the analysis.
384   AliInfo( Form( "Processed %i tracks.", fTrackTot ) );
385   AliInfo( Form( "Accepted %i tracks.", fTrackOk ) );
386
387   // stop here if no alignment is to be performed
388   if( !fDoAlignment ) return;
389
390   AliLog::SetGlobalLogLevel(AliLog::kInfo);
391
392   // Perform global fit
393   fAlign->GlobalFit(fParameters,fErrors,fPulls);
394   Int_t idOffset = 0; // 400
395   Int_t lSDetElemCh = 0;
396
397   for( Int_t iDE=0; iDE<AliMUONAlignment::fgNDetElem; iDE++ )
398   {
399
400     Int_t detElemId = idOffset+100;
401     detElemId += iDE;
402     lSDetElemCh = 0;
403     for(Int_t iCh=0; iCh<9; iCh++)
404     {
405
406       lSDetElemCh += AliMUONAlignment::fgNDetElemCh[iCh];
407       if (iDE>=lSDetElemCh)
408       {
409         detElemId += 100;
410         detElemId -= AliMUONAlignment::fgNDetElemCh[iCh];
411       }
412
413     }
414
415   }
416
417   // Post final data. Write histo list to a file with option "RECREATE"
418   // PostData(1,fList);
419
420   // store misalignments from OCDB into old transformers
421   if( fMergeAlignmentCDBs )
422   { SaveMisAlignmentData( fOldGeoTransformer ); }
423
424   // Re Align
425   fNewGeoTransformer = fAlign->ReAlign( fOldGeoTransformer, fParameters, true );
426
427   // Generate realigned data in local cdb
428   const TClonesArray* array = fNewGeoTransformer->GetMisAlignmentData();
429
430   // 100 mum residual resolution for chamber misalignments?
431   fAlign->SetAlignmentResolution( array, -1, 0.01, 0.01, 0.004, 0.003 );
432
433   // CDB manager
434   AliLog::SetGlobalDebugLevel(2);
435
436   // recover default storage full name (raw:// cannot be used to set specific storage)
437   AliCDBManager* cdbManager = AliCDBManager::Instance();
438
439   // unload old alignment path
440   if( cdbManager->GetEntryCache()->Contains("MUON/Align/Data") )
441   { cdbManager->UnloadFromCache("MUON/Align/Data"); }
442
443   // load new alignment path
444   if( !fNewAlignStorage.IsNull() ) cdbManager->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
445   else {
446
447     TString defaultStorage(cdbManager->GetDefaultStorage()->GetType());
448     if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
449     else defaultStorage += Form("://%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
450     cdbManager->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
451
452   }
453
454   // create new DB entry
455   AliCDBMetaData* cdbData = new AliCDBMetaData();
456   cdbData->SetResponsible("Dimuon Offline project");
457   cdbData->SetComment("MUON alignment objects with residual misalignment");
458   AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
459   cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
460
461 }
462
463 //________________________________________________________________________
464 void AliMUONAlignmentTask::NotifyRun()
465 {
466
467   /// run number (re)initialization
468
469   // propagate run number to fAlign
470   if( fAlign ) fAlign->SetRunNumber( fCurrentRunNumber );
471
472   // update ocdb
473   if (fOCDBLoaded && fLoadOCDBOnce)
474   {
475     AliError(Form("OCDB already loaded %d %d",fOCDBLoaded,fLoadOCDBOnce));
476     return;
477   }
478
479   AliCDBManager* cdbManager = AliCDBManager::Instance();
480
481   // Initialize default storage
482   if( !( cdbManager->IsDefaultStorageSet() || fDefaultStorage.IsNull() ) )
483   {
484
485     AliInfo( Form( "Initializing default storage: %s", fDefaultStorage.Data() ) );
486     cdbManager->SetDefaultStorage(fDefaultStorage.Data());
487
488   } else if( !fDefaultStorage.IsNull() ) {
489
490     AliInfo( "Default storage already set. Ignoring fDefaultStorage" );
491
492   } else {
493
494     AliInfo( "Default storage already set" );
495
496   }
497
498   // Initialize run number
499   if( cdbManager->GetRun() < 0 )
500   {
501     AliInfo( Form( "Setting run number: %i", fCurrentRunNumber ) );
502     cdbManager->SetRun(fCurrentRunNumber);
503   }
504
505   // following initialization is not needed when reading records
506   if( !fReadRecords )
507   {
508
509     // load magnetic field if needed
510     if( !TGeoGlobalMagField::Instance()->IsLocked() )
511     {
512
513       AliInfo( "Loading magnetic field" );
514       if( !AliMUONCDB::LoadField() )
515       {
516         AliError( "Failed to load magnetic field" );
517         return;
518       }
519
520     } else { AliInfo( "Magnetic field already locked" ); }
521
522     // checking magnitic field
523     if( !fForceBField )
524     {
525       AliInfo( "Reading magnetic field setup from GRP" );
526
527       // decide bFieldOn value base on dipole current
528       // propagete to Alignment class
529       // and printout
530       AliMagF* magF = dynamic_cast<AliMagF*>( TGeoGlobalMagField::Instance()->GetField() );
531       fBFieldOn = TMath::Abs( magF->GetFactorDip() ) > 1e-5;
532       fAlign->SetBFieldOn( fBFieldOn );
533       AliInfo( Form( "Dipole magnetic field factor: %.2f", magF->GetFactorDip() ) );
534       AliInfo( Form( "fBFieldOn = %s", (fBFieldOn ? "kTRUE":"kFALSE") ) );
535     }
536
537     AliInfo( "Loading muon mapping" );
538     if( !AliMUONCDB::LoadMapping(kFALSE) )
539     {
540       AliError( "Failed to load muon mapping" );
541       return;
542     }
543
544     AliInfo( "Assigning field to Track extrapolator" );
545     AliMUONTrackExtrap::SetField();
546
547   }
548
549   // load geometry if needed
550   if( !AliGeomManager::GetGeometry() )
551   {
552
553     // reset existing geometry/alignment if any
554     if( cdbManager->GetEntryCache()->Contains("GRP/Geometry/Data") )
555     {
556       AliInfo( "Unloading GRP/Geometry/Data" );
557       cdbManager->UnloadFromCache("GRP/Geometry/Data");
558     }
559
560     if( cdbManager->GetEntryCache()->Contains("MUON/Align/Data") )
561     {
562       AliInfo( Form( "Unloading MUON/Align/Data from %s", cdbManager->GetSpecificStorage( "MUON/Align/Data" )->GetBaseFolder().Data() ) );
563       cdbManager->UnloadFromCache("MUON/Align/Data");
564     }
565
566     // get original geometry transformer
567     AliInfo( "Loading geometry" );
568     AliGeomManager::LoadGeometry();
569     if (!AliGeomManager::GetGeometry())
570     {
571       AliError("Failed to load geometry");
572       return;
573     }
574
575     if (!fOldAlignStorage.IsNull())
576     {
577
578       AliInfo( Form( "Initializing MUON/Align/Data using: %s", fOldAlignStorage.Data() ) );
579       cdbManager->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
580
581     } else {
582
583       // recover default storage full name (raw:// cannot be used to set specific storage)
584       TString defaultStorage(cdbManager->GetDefaultStorage()->GetType());
585       if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
586       else defaultStorage += Form("://%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
587
588       AliInfo( Form( "Re-initializing MUON/Align/Data using: %s", defaultStorage.Data() ) );
589       cdbManager->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
590
591     }
592
593     AliInfo("Loading muon Alignment objects");
594     AliGeomManager::ApplyAlignObjsFromCDB("MUON");
595
596   } else { AliInfo( "Geometry already initialized - fOldAlignStorage ignored" ); }
597
598   // update geometry transformer and pass to Alignment object
599   AliInfo("Loading muon geometry in transformer");
600   fOldGeoTransformer->LoadGeometryData();
601   fAlign->SetGeometryTransformer(fOldGeoTransformer);
602
603   fOCDBLoaded = kTRUE;
604
605 }
606
607 //_____________________________________________________________________________________
608 void AliMUONAlignmentTask::SaveMisAlignmentData( AliMUONGeometryTransformer* transformer ) const
609 {
610
611   // clear transformer
612   transformer->ClearMisAlignmentData();
613
614   // load MUON/Align/Data from OCDB
615   AliCDBManager* cdbManager = AliCDBManager::Instance();
616   AliCDBEntry* cdbEntry = cdbManager->Get( "MUON/Align/Data" );
617   if (!cdbEntry)
618   {
619
620     AliError( "unable to load entry for path MUON/Align/Data" );
621     return;
622
623   }
624
625   // get TClonesArray and check
626   TClonesArray* misArray = (TClonesArray*)cdbEntry->GetObject();
627   if( !misArray )
628   {
629
630     AliError( "unable to load old misalignment array" );
631     return;
632
633   }
634
635   // loop over stored entries
636   for (Int_t index=0; index<misArray->GetEntriesFast(); ++index )
637   {
638
639     // load matrix and check
640     AliAlignObjMatrix* matrix = (AliAlignObjMatrix*) misArray->At( index );
641     if( !matrix )
642     {
643       AliError( Form( "unable to load matrix for index %i", index ) );
644       continue;
645     }
646
647     // get volume ID
648     AliGeomManager::ELayerID layerId;
649     Int_t moduleId;
650     matrix->GetVolUID( layerId, moduleId);
651
652     // make sure ELayerID is correct
653     assert( layerId == AliGeomManager::kMUON );
654
655     // printout
656     // AliInfo( Form( "Found matrix for %s %i", matrix->GetSymName(), moduleId ) );
657
658     // get matrix
659     TGeoHMatrix misMatrix;
660     matrix->GetMatrix(misMatrix);
661
662     // add to geometry transformer
663     // need the detector element
664     // "detElement->GetId()"
665     // see fOldGeoTransformer->GetMisAlignmentData( ... )
666
667     if( TString( matrix->GetSymName() ).Contains( "DE" ) ) transformer->AddMisAlignDetElement( moduleId,  misMatrix);
668     else transformer->AddMisAlignModule( moduleId,  misMatrix);
669   }
670
671   return;
672 }
673