Fixes to run the MUON alignment on the grid
[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 AliMillepede.
23 /// A OCDB entry is written with the alignment parameters.
24 /// A list of graph are written to monitor the alignment parameters.
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 <TGraphErrors.h>
34 #include <TTree.h>
35 #include <TChain.h>
36 #include <TClonesArray.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 "AliCDBMetaData.h"
50 #include "AliCDBId.h"
51 #include "AliGeomManager.h"
52
53 #include "AliMpCDB.h"
54 #include "AliMUONAlignment.h"
55 #include "AliMUONTrack.h"
56 #include "AliMUONTrackExtrap.h"
57 #include "AliMUONTrackParam.h"
58 #include "AliMUONGeometryTransformer.h"
59 #include "AliMUONESDInterface.h"
60
61 #include "AliMUONAlignmentTask.h"
62
63 ///\cond CLASSIMP  
64 ClassImp(AliMUONAlignmentTask)
65 ///\endcond
66
67 // //________________________________________________________________________
68 // AliMUONAlignmentTask::AliMUONAlignmentTask(const char *name) 
69 //   : AliAnalysisTask(name, ""),
70 //     fESD(0x0),
71 //     fAlign(0x0),
72 //     fGeoFilename(""),
73 //     fTransform(0x0),
74 //     fTrackTot(0),
75 //     fTrackOk(0),
76 //     fMSDEx(0x0), 
77 //     fMSDEy(0x0), 
78 //     fMSDEz(0x0),
79 //     fMSDEp(0x0),
80 //     fList(0x0)
81 // {
82 //   /// Default Constructor
83 //   // Define input and output slots here
84 //   // Input slot #0 works with a TChain
85 //   DefineInput(0, TChain::Class());
86 //   // Output slot #0 writes NTuple/histos into a TList
87 //   DefineOutput(0, TList::Class());  
88
89 //   // initialize parameters ...
90 //   for(Int_t k=0;k<4*156;k++) {
91 //     fParameters[k]=0.;
92 //     fErrors[k]=0.;
93 //     fPulls[k]=0.;
94 //   }
95
96 //   fAlign = new AliMUONAlignment();
97 //   fTransform = new AliMUONGeometryTransformer();  
98 // }
99
100 //________________________________________________________________________
101 AliMUONAlignmentTask::AliMUONAlignmentTask(const char *name, const char *geofilename, const char *defaultocdb, const char *misalignocdb) 
102   : AliAnalysisTask(name, ""),
103     fESD(0x0),
104     fAlign(0x0),
105     fGeoFilename(geofilename),
106     fMisAlignOCDB(misalignocdb),
107     fDefaultOCDB(defaultocdb),
108     fTransform(0x0),
109     fTrackTot(0),
110     fTrackOk(0),
111     fLastRunNumber(-1),
112     fMSDEx(0x0), 
113     fMSDEy(0x0), 
114     fMSDEz(0x0),
115     fMSDEp(0x0),
116     fList(0x0)
117 {
118   /// Default Constructor
119   // Define input and output slots here
120   // Input slot #0 works with a TChain
121   DefineInput(0, TChain::Class());
122   // Output slot #0 writes NTuple/histos into a TList
123   DefineOutput(0, TList::Class());  
124
125   // initialize parameters ...
126   for(Int_t k=0;k<4*156;k++) {
127     fParameters[k]=0.;
128     fErrors[k]=0.;
129     fPulls[k]=0.;
130   }
131
132   fAlign = new AliMUONAlignment();
133   fTransform = new AliMUONGeometryTransformer();  
134 }
135
136
137 //________________________________________________________________________
138 AliMUONAlignmentTask::AliMUONAlignmentTask(const AliMUONAlignmentTask& obj) 
139   : AliAnalysisTask(obj),
140     fESD(0x0),
141     fAlign(0x0),
142     fGeoFilename(""),
143     fMisAlignOCDB(""),
144     fDefaultOCDB(""),
145     fTransform(0x0),
146     fTrackTot(0),
147     fTrackOk(0),
148     fLastRunNumber(-1),
149     fMSDEx(0x0), 
150     fMSDEy(0x0), 
151     fMSDEz(0x0),
152     fMSDEp(0x0),
153     fList(0x0)
154 {
155   /// Copy constructor
156   fESD = obj.fESD;
157   fAlign = obj.fAlign;
158   fGeoFilename = obj.fGeoFilename;
159   fTransform = obj.fTransform;
160   fTrackTot = obj.fTrackTot;  
161   fTrackOk = obj.fTrackOk;  
162   fLastRunNumber = obj.fLastRunNumber;
163   fMSDEx = obj.fMSDEx; 
164   fMSDEy = obj.fMSDEy; 
165   fMSDEz = obj.fMSDEz;
166   fMSDEp = obj.fMSDEp;
167   fList = obj.fList;  
168 }
169
170 //________________________________________________________________________
171 AliMUONAlignmentTask& AliMUONAlignmentTask::operator=(const AliMUONAlignmentTask& other) 
172 {
173   /// Assignment
174   AliAnalysisTask::operator=(other);
175   fESD = other.fESD;
176   fAlign = other.fAlign;
177   fGeoFilename = other.fGeoFilename;
178   fMisAlignOCDB = other.fMisAlignOCDB;
179   fDefaultOCDB = other.fDefaultOCDB;
180   fTransform = other.fTransform;
181   fTrackTot = other.fTrackTot;  
182   fTrackOk = other.fTrackOk;  
183   fLastRunNumber = other.fLastRunNumber;
184   fMSDEx = other.fMSDEx; 
185   fMSDEy = other.fMSDEy; 
186   fMSDEz = other.fMSDEz;
187   fMSDEp = other.fMSDEp;
188   fList = other.fList;  
189   
190   return *this;
191 }
192
193 //________________________________________________________________________
194 AliMUONAlignmentTask::~AliMUONAlignmentTask() 
195
196   /// Destructor
197   if (fAlign) delete fAlign;
198   if (fTransform) delete fTransform;
199 }
200
201 //________________________________________________________________________
202 void AliMUONAlignmentTask::LocalInit() 
203 {
204   /// Local initialization, called once per task on the client machine 
205   /// where the analysis train is assembled
206   fLastRunNumber = 0; 
207   //  Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
208   Prepare(fGeoFilename.Data(),"local://$ALICE_ROOT/OCDB",fMisAlignOCDB.Data());
209   fLastRunNumber = -1;
210
211   // Set initial values here, good guess may help convergence
212   // St 1 
213   //  Int_t iPar = 0;
214   //  fParameters[iPar++] =  0.010300 ;  fParameters[iPar++] =  0.010600 ;  fParameters[iPar++] =  0.000396 ;  
215
216
217   fAlign->InitGlobalParameters(fParameters);
218
219
220   fTransform->LoadGeometryData();
221   fAlign->SetGeometryTransformer(fTransform);
222
223   // Do alignment with magnetic field off
224   fAlign->SetBFieldOn(kFALSE);
225   
226   // Set tracking station to use
227   //  Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
228   Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kFALSE,kTRUE,kTRUE,kTRUE,kTRUE};
229
230   // Set degrees of freedom to align (see AliMUONAlignment)
231   fAlign->AllowVariations(bChOnOff);
232
233   // Fix parameters or add constraints here
234   //   for (Int_t iSt=0; iSt<5; iSt++)
235   //     if (!bStOnOff[iSt]) fAlign->FixStation(iSt+1);
236   for (Int_t iCh=0; iCh<10; iCh++)
237     if (!bChOnOff[iCh]) fAlign->FixChamber(iCh+1);
238
239   // Left and right sides of the detector are independent, one can choose to align 
240   // only one side
241   Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
242   fAlign->FixHalfSpectrometer(bChOnOff,bSpecLROnOff);
243
244   fAlign->SetChOnOff(bChOnOff);
245   fAlign->SetSpecLROnOff(bChOnOff);
246
247     // Here we can fix some detection elements
248   fAlign->FixDetElem(908);
249   fAlign->FixDetElem(1020);
250
251   // Set predifined global constrains: X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
252 //   Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
253 //   Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
254   //  fAlign->AddConstraints(bChOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
255
256 }
257
258 //________________________________________________________________________
259 void AliMUONAlignmentTask::ConnectInputData(Option_t *) 
260 {
261   /// Connect ESD here. Called on each input data change.
262
263   // Connect ESD here
264   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
265   if (!esdTree) {
266     Printf("ERROR: Could not read chain from input slot 0");
267   } 
268   else {
269     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());   
270     if (!esdH) {
271       Printf("ERROR: Could not get ESDInputHandler");
272     } 
273     else {
274       fESD = esdH->GetEvent();
275     }
276   }
277 }
278
279 //________________________________________________________________________
280 void AliMUONAlignmentTask::CreateOutputObjects()
281 {
282   /// Executed once on each worker (machine actually running the analysis code)
283   //
284   // This method has to be called INSIDE the user redefined CreateOutputObjects
285   // method, before creating each object corresponding to the output containers
286   // that are to be written to a file. This need to be done in general for the big output
287   // objects that may not fit memory during processing. 
288   //  OpenFile(0); 
289
290   // Creating graphs 
291   fMSDEx = new TGraphErrors(156);
292   fMSDEy = new TGraphErrors(156);
293   fMSDEz = new TGraphErrors(156);
294   fMSDEp = new TGraphErrors(156);
295
296   // Add Ntuples to the list
297   fList = new TList();
298   fList->Add(fMSDEx);
299   fList->Add(fMSDEy);
300   fList->Add(fMSDEz);
301   fList->Add(fMSDEp);
302 }
303
304 //________________________________________________________________________
305 void AliMUONAlignmentTask::Exec(Option_t *) 
306 {
307   /// Main loop, called for each event
308   if (!fESD) {
309     Printf("ERROR: fESD not available");
310     return;
311   }
312
313   Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
314   printf("ESD RunNumber %d , fLastRunNumber %d\n",fESD->GetRunNumber(),fLastRunNumber);
315   if (fESD->GetRunNumber()!=fLastRunNumber){
316     fLastRunNumber = fESD->GetRunNumber();
317     Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
318   }
319  
320   Int_t nTracks = Int_t(fESD->GetNumberOfMuonTracks());
321   //  if (!event%100) cout << " there are " << nTracks << " tracks in event " << event << endl;
322   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
323     AliESDMuonTrack* esdTrack = fESD->GetMuonTrack(iTrack);
324     if (!esdTrack->ClustersStored()) continue;
325     Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
326 //     fInvBenMom->Fill(invBenMom);
327 //     fBenMom->Fill(1./invBenMom);
328     if (TMath::Abs(invBenMom)<=1.04) {
329       AliMUONTrack track;
330       printf("cero\n");
331       AliMUONESDInterface::ESDToMUON(*esdTrack, track);
332       printf("uno %p\n", fAlign);
333       fAlign->ProcessTrack(&track);
334       printf("dos\n");
335       fAlign->LocalFit(fTrackOk++,trackParams,0);
336       printf("tres\n");
337     }
338     fTrackTot++;
339     cout << "Processed " << fTrackTot << " Tracks." << endl;
340   }
341   
342   // Post final data. Write histo list to a file with option "RECREATE"
343   PostData(0,fList);
344 }      
345
346 //________________________________________________________________________
347 void AliMUONAlignmentTask::FinishTaskOutput()
348 {
349   /// Called once per task on the client machine at the end of the analysis.
350
351   cout << "Processed " << fTrackTot << " Tracks." << endl;
352   // Perform global fit
353   fAlign->GlobalFit(fParameters,fErrors,fPulls);
354
355   cout << "Done with GlobalFit " << endl;
356
357 //   // Update pointers reading them from the output slot
358 //   fList = (TList*)GetOutputData(0);
359 //   fMSDEx = (TGraphErrors*)fList->At(0);
360 //   fMSDEy = (TGraphErrors*)fList->At(1);
361 //   fMSDEz = (TGraphErrors*)fList->At(2);
362 //   fMSDEp = (TGraphErrors*)fList->At(3);
363
364   // Store results
365   Double_t DEid[156] = {0};
366   Double_t MSDEx[156] = {0};
367   Double_t MSDEy[156] = {0};
368   Double_t MSDEz[156] = {0};
369   Double_t MSDEp[156] = {0};
370   Double_t DEidErr[156] = {0};
371   Double_t MSDExErr[156] = {0};
372   Double_t MSDEyErr[156] = {0};
373   Double_t MSDEzErr[156] = {0};
374   Double_t MSDEpErr[156] = {0};
375   Int_t lNDetElem = 4*2+4*2+18*2+26*2+26*2;
376   Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
377   // Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
378   Int_t idOffset = 0; // 400
379   Int_t lSDetElemCh = 0;
380   for(Int_t iDE=0; iDE<lNDetElem; iDE++){
381     DEidErr[iDE] = 0.;
382     DEid[iDE] = idOffset+100;
383     DEid[iDE] += iDE; 
384     lSDetElemCh = 0;
385     for(Int_t iCh=0; iCh<9; iCh++){
386       lSDetElemCh += lNDetElemCh[iCh];
387       if (iDE>=lSDetElemCh) {
388         DEid[iDE] += 100;
389         DEid[iDE] -= lNDetElemCh[iCh];
390       }
391     }
392     MSDEx[iDE]=fParameters[3*iDE+0];
393     MSDEy[iDE]=fParameters[3*iDE+1];
394     MSDEz[iDE]=fParameters[3*iDE+3];
395     MSDEp[iDE]=fParameters[3*iDE+2];
396     MSDExErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+0);
397     MSDEyErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+1);
398     MSDEzErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+3);
399     MSDEpErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+2);
400     fMSDEx->SetPoint(iDE,DEid[iDE],fParameters[3*iDE+0]);
401     fMSDEx->SetPoint(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(3*iDE+0));
402     fMSDEy->SetPoint(iDE,DEid[iDE],fParameters[3*iDE+1]);
403     fMSDEy->SetPoint(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(3*iDE+1));
404     fMSDEp->SetPoint(iDE,DEid[iDE],fParameters[3*iDE+2]);
405     fMSDEp->SetPoint(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(3*iDE+2));
406     fMSDEz->SetPoint(iDE,DEid[iDE],fParameters[3*iDE+3]);
407     fMSDEz->SetPoint(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(3*iDE+3));
408   }
409
410   // Post final data. Write histo list to a file with option "RECREATE"
411   PostData(0,fList);
412
413   // Re Align
414   AliMUONGeometryTransformer *newTransform = fAlign->ReAlign(fTransform,fParameters,true); 
415   newTransform->WriteTransformations("transform2ReAlign.dat");
416   
417   // Generate realigned data in local cdb
418   const TClonesArray* array = newTransform->GetMisAlignmentData();
419
420   // 100 mum residual resolution for chamber misalignments?
421   fAlign->SetAlignmentResolution(array,-1,0.01,0.01,0.004,0.003);
422    
423   // CDB manager
424   AliCDBManager* cdbManager = AliCDBManager::Instance();
425   cdbManager->SetDefaultStorage(fDefaultOCDB.Data());
426   cdbManager->SetSpecificStorage("MUON/Align/Data",fMisAlignOCDB.Data());
427   
428   AliCDBMetaData* cdbData = new AliCDBMetaData();
429   cdbData->SetResponsible("Dimuon Offline project");
430   cdbData->SetComment("MUON alignment objects with residual misalignment");
431   AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity()); 
432   cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
433
434 }
435
436 //________________________________________________________________________
437 void AliMUONAlignmentTask::Terminate(const Option_t*)
438 {
439   /// Called once per task on the client machine at the end of the analysis.
440
441 }
442
443 //-----------------------------------------------------------------------
444 void AliMUONAlignmentTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
445 {
446   /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
447       
448   // Load mapping
449   AliCDBManager* man = AliCDBManager::Instance();
450   man->SetDefaultStorage(defaultOCDB);
451   man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
452   man->Print();
453   man->SetRun(fLastRunNumber);
454   if ( ! AliMpCDB::LoadDDLStore() ) {
455     Error("MUONRefit","Could not access mapping from OCDB !");
456     exit(-1);
457   }
458
459   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
460   if (!gGeoManager) {
461     AliGeomManager::LoadGeometry(geoFilename);
462     if (!gGeoManager) {
463       Error("AliMUONReAlignTask", "getting geometry from file %s failed", "generated/galice.root");
464       return;
465     }
466   }
467
468   // set mag field
469   if (!TGeoGlobalMagField::Instance()->GetField()) {
470     printf("Loading field map...\n");
471     AliGRPManager *grpMan = new AliGRPManager();
472     grpMan->ReadGRPEntry();
473     grpMan->SetMagField();
474     delete grpMan;
475   }
476   // set the magnetic field for track extrapolations
477   AliMUONTrackExtrap::SetField();
478   
479 }