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