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