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   if (fESD->GetRunNumber()!=fLastRunNumber){
 315     fLastRunNumber = fESD->GetRunNumber();
 316     Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
 317   }
 318
 319   Int_t nTracks = Int_t(fESD->GetNumberOfMuonTracks());
 320   //  if (!event%100) cout << " there are " << nTracks << " tracks in event " << event << endl;
 321   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 322     AliESDMuonTrack* esdTrack = fESD->GetMuonTrack(iTrack);
 323     if (!esdTrack->ClustersStored()) continue;
 324     Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
 325 //     fInvBenMom->Fill(invBenMom);
 326 //     fBenMom->Fill(1./invBenMom);
 327     if (TMath::Abs(invBenMom)<=1.04) {
 328       AliMUONTrack track;
 329       AliMUONESDInterface::ESDToMUON(*esdTrack, track);
 330       fAlign->ProcessTrack(&track);
 331       fAlign->LocalFit(fTrackOk++,trackParams,0);
 332     }
 333     fTrackTot++;
 334     cout << "Processed " << fTrackTot << " Tracks." << endl;
 335   }
 336
 337   // Post final data. Write histo list to a file with option "RECREATE"
 338   PostData(0,fList);
 339 }
 340
 341 //________________________________________________________________________
 342 void AliMUONAlignmentTask::FinishTaskOutput()
 343 {
 344   /// Called once per task on the client machine at the end of the analysis.
 345
 346   cout << "Processed " << fTrackTot << " Tracks." << endl;
 347   // Perform global fit
 348   fAlign->GlobalFit(fParameters,fErrors,fPulls);
 349
 350   cout << "Done with GlobalFit " << endl;
 351
 352 //   // Update pointers reading them from the output slot
 353 //   fList = (TList*)GetOutputData(0);
 354 //   fMSDEx = (TGraphErrors*)fList->At(0);
 355 //   fMSDEy = (TGraphErrors*)fList->At(1);
 356 //   fMSDEz = (TGraphErrors*)fList->At(2);
 357 //   fMSDEp = (TGraphErrors*)fList->At(3);
 358
 359   // Store results
 360   Double_t deId[156] = {0};
 361   Double_t msdEx[156] = {0};
 362   Double_t msdEy[156] = {0};
 363   Double_t msdEz[156] = {0};
 364   Double_t msdEp[156] = {0};
 365   Double_t deIdErr[156] = {0};
 366   Double_t msdExErr[156] = {0};
 367   Double_t msdEyErr[156] = {0};
 368   Double_t msdEzErr[156] = {0};
 369   Double_t msdEpErr[156] = {0};
 370   Int_t lNDetElem = 4*2+4*2+18*2+26*2+26*2;
 371   Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
 372   // Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
 373   Int_t idOffset = 0; // 400
 374   Int_t lSDetElemCh = 0;
 375   for(Int_t iDE=0; iDE<lNDetElem; iDE++){
 376     deIdErr[iDE] = 0.;
 377     deId[iDE] = idOffset+100;
 378     deId[iDE] += iDE;
 379     lSDetElemCh = 0;
 380     for(Int_t iCh=0; iCh<9; iCh++){
 381       lSDetElemCh += lNDetElemCh[iCh];
 382       if (iDE>=lSDetElemCh) {
 383         deId[iDE] += 100;
 384         deId[iDE] -= lNDetElemCh[iCh];
 385       }
 386     }
 387     msdEx[iDE]=fParameters[3*iDE+0];
 388     msdEy[iDE]=fParameters[3*iDE+1];
 389     msdEz[iDE]=fParameters[3*iDE+3];
 390     msdEp[iDE]=fParameters[3*iDE+2];
 391     msdExErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+0);
 392     msdEyErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+1);
 393     msdEzErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+3);
 394     msdEpErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+2);
 395     fMSDEx->SetPoint(iDE,deId[iDE],fParameters[3*iDE+0]);
 396     fMSDEx->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+0));
 397     fMSDEy->SetPoint(iDE,deId[iDE],fParameters[3*iDE+1]);
 398     fMSDEy->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+1));
 399     fMSDEp->SetPoint(iDE,deId[iDE],fParameters[3*iDE+2]);
 400     fMSDEp->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+2));
 401     fMSDEz->SetPoint(iDE,deId[iDE],fParameters[3*iDE+3]);
 402     fMSDEz->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+3));
 403   }
 404
 405   // Post final data. Write histo list to a file with option "RECREATE"
 406   PostData(0,fList);
 407
 408   // Re Align
 409   AliMUONGeometryTransformer *newTransform = fAlign->ReAlign(fTransform,fParameters,true);
 410   newTransform->WriteTransformations("transform2ReAlign.dat");
 411
 412   // Generate realigned data in local cdb
 413   const TClonesArray* array = newTransform->GetMisAlignmentData();
 414
 415   // 100 mum residual resolution for chamber misalignments?
 416   fAlign->SetAlignmentResolution(array,-1,0.01,0.01,0.004,0.003);
 417
 418   // CDB manager
 419   AliCDBManager* cdbManager = AliCDBManager::Instance();
 420   cdbManager->SetDefaultStorage(fDefaultOCDB.Data());
 421   cdbManager->SetSpecificStorage("MUON/Align/Data",fMisAlignOCDB.Data());
 422
 423   AliCDBMetaData* cdbData = new AliCDBMetaData();
 424   cdbData->SetResponsible("Dimuon Offline project");
 425   cdbData->SetComment("MUON alignment objects with residual misalignment");
 426   AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
 427   cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
 428
 429 }
 430
 431 //________________________________________________________________________
 432 void AliMUONAlignmentTask::Terminate(const Option_t*)
 433 {
 434   /// Called once per task on the client machine at the end of the analysis.
 435
 436 }
 437
 438 //-----------------------------------------------------------------------
 439 void AliMUONAlignmentTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
 440 {
 441   /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
 442
 443   // Load mapping
 444   AliCDBManager* man = AliCDBManager::Instance();
 445   man->SetDefaultStorage(defaultOCDB);
 446   man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
 447   man->Print();
 448   man->SetRun(fLastRunNumber);
 449   if ( ! AliMpCDB::LoadDDLStore() ) {
 450     Error("MUONRefit","Could not access mapping from OCDB !");
 451     exit(-1);
 452   }
 453
 454   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
 455   if (!gGeoManager) {
 456     AliGeomManager::LoadGeometry(geoFilename);
 457     if (!gGeoManager) {
 458       Error("AliMUONReAlignTask", "getting geometry from file %s failed", "generated/galice.root");
 459       return;
 460     }
 461   }
 462
 463   // set mag field
 464   if (!TGeoGlobalMagField::Instance()->GetField()) {
 465     printf("Loading field map...\n");
 466     AliGRPManager *grpMan = new AliGRPManager();
 467     grpMan->ReadGRPEntry();
 468     grpMan->SetMagField();
 469     delete grpMan;
 470   }
 471   // set the magnetic field for track extrapolations
 472   AliMUONTrackExtrap::SetField();
 473
 474 }