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   // initialize parameters ...
 170   for(Int_t k=0;k<4*156;k++) {
 171     fParameters[k]=obj.fParameters[k];
 172     fErrors[k]=obj.fErrors[k];
 173     fPulls[k]=obj.fPulls[k];
 174   }
 175
 176 }
 177
 178 //________________________________________________________________________
 179 AliMUONAlignmentTask& AliMUONAlignmentTask::operator=(const AliMUONAlignmentTask& other)
 180 {
 181   /// Assignment
 182   AliAnalysisTask::operator=(other);
 183   fESD = other.fESD;
 184   fAlign = other.fAlign;
 185   fGeoFilename = other.fGeoFilename;
 186   fMisAlignOCDB = other.fMisAlignOCDB;
 187   fDefaultOCDB = other.fDefaultOCDB;
 188   fTransform = other.fTransform;
 189   fTrackTot = other.fTrackTot;
 190   fTrackOk = other.fTrackOk;
 191   fLastRunNumber = other.fLastRunNumber;
 192   fMSDEx = other.fMSDEx;
 193   fMSDEy = other.fMSDEy;
 194   fMSDEz = other.fMSDEz;
 195   fMSDEp = other.fMSDEp;
 196   fList = other.fList;
 197
 198   // initialize parameters ...
 199   for(Int_t k=0;k<4*156;k++) {
 200     fParameters[k]=other.fParameters[k];
 201     fErrors[k]=other.fErrors[k];
 202     fPulls[k]=other.fPulls[k];
 203   }
 204
 205   return *this;
 206 }
 207
 208 //________________________________________________________________________
 209 AliMUONAlignmentTask::~AliMUONAlignmentTask()
 210 {
 211   /// Destructor
 212   if (fAlign) delete fAlign;
 213   if (fTransform) delete fTransform;
 214 }
 215
 216 //________________________________________________________________________
 217 void AliMUONAlignmentTask::LocalInit()
 218 {
 219   /// Local initialization, called once per task on the client machine
 220   /// where the analysis train is assembled
 221   fLastRunNumber = 0;
 222   //  Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
 223   Prepare(fGeoFilename.Data(),"local://$ALICE_ROOT/OCDB",fMisAlignOCDB.Data());
 224   fLastRunNumber = -1;
 225
 226   // Set initial values here, good guess may help convergence
 227   // St 1
 228   //  Int_t iPar = 0;
 229   //  fParameters[iPar++] =  0.010300 ;  fParameters[iPar++] =  0.010600 ;  fParameters[iPar++] =  0.000396 ;
 230
 231
 232   fAlign->InitGlobalParameters(fParameters);
 233
 234
 235   fTransform->LoadGeometryData();
 236   fAlign->SetGeometryTransformer(fTransform);
 237
 238   // Do alignment with magnetic field off
 239   fAlign->SetBFieldOn(kFALSE);
 240
 241   // Set tracking station to use
 242   //  Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
 243   Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kFALSE,kTRUE,kTRUE,kTRUE,kTRUE};
 244
 245   // Set degrees of freedom to align (see AliMUONAlignment)
 246   fAlign->AllowVariations(bChOnOff);
 247
 248   // Fix parameters or add constraints here
 249   //   for (Int_t iSt=0; iSt<5; iSt++)
 250   //     if (!bStOnOff[iSt]) fAlign->FixStation(iSt+1);
 251   for (Int_t iCh=0; iCh<10; iCh++)
 252     if (!bChOnOff[iCh]) fAlign->FixChamber(iCh+1);
 253
 254   // Left and right sides of the detector are independent, one can choose to align
 255   // only one side
 256   Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
 257   fAlign->FixHalfSpectrometer(bChOnOff,bSpecLROnOff);
 258
 259   fAlign->SetChOnOff(bChOnOff);
 260   fAlign->SetSpecLROnOff(bChOnOff);
 261
 262     // Here we can fix some detection elements
 263   fAlign->FixDetElem(908);
 264   fAlign->FixDetElem(1020);
 265
 266   // Set predifined global constrains: X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
 267 //   Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
 268 //   Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
 269   //  fAlign->AddConstraints(bChOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
 270
 271 }
 272
 273 //________________________________________________________________________
 274 void AliMUONAlignmentTask::ConnectInputData(Option_t *)
 275 {
 276   /// Connect ESD here. Called on each input data change.
 277
 278   // Connect ESD here
 279   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
 280   if (!esdTree) {
 281     Printf("ERROR: Could not read chain from input slot 0");
 282   }
 283   else {
 284     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 285     if (!esdH) {
 286       Printf("ERROR: Could not get ESDInputHandler");
 287     }
 288     else {
 289       fESD = esdH->GetEvent();
 290     }
 291   }
 292 }
 293
 294 //________________________________________________________________________
 295 void AliMUONAlignmentTask::CreateOutputObjects()
 296 {
 297   /// Executed once on each worker (machine actually running the analysis code)
 298   //
 299   // This method has to be called INSIDE the user redefined CreateOutputObjects
 300   // method, before creating each object corresponding to the output containers
 301   // that are to be written to a file. This need to be done in general for the big output
 302   // objects that may not fit memory during processing.
 303   //  OpenFile(0);
 304
 305   // Creating graphs
 306   fMSDEx = new TGraphErrors(156);
 307   fMSDEy = new TGraphErrors(156);
 308   fMSDEz = new TGraphErrors(156);
 309   fMSDEp = new TGraphErrors(156);
 310
 311   // Add Ntuples to the list
 312   fList = new TList();
 313   fList->Add(fMSDEx);
 314   fList->Add(fMSDEy);
 315   fList->Add(fMSDEz);
 316   fList->Add(fMSDEp);
 317 }
 318
 319 //________________________________________________________________________
 320 void AliMUONAlignmentTask::Exec(Option_t *)
 321 {
 322   /// Main loop, called for each event
 323   if (!fESD) {
 324     Printf("ERROR: fESD not available");
 325     return;
 326   }
 327
 328   Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
 329   if (fESD->GetRunNumber()!=fLastRunNumber){
 330     fLastRunNumber = fESD->GetRunNumber();
 331     Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
 332   }
 333
 334   Int_t nTracks = Int_t(fESD->GetNumberOfMuonTracks());
 335   //  if (!event%100) cout << " there are " << nTracks << " tracks in event " << event << endl;
 336   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 337     AliESDMuonTrack* esdTrack = fESD->GetMuonTrack(iTrack);
 338     if (!esdTrack->ClustersStored()) continue;
 339     Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
 340 //     fInvBenMom->Fill(invBenMom);
 341 //     fBenMom->Fill(1./invBenMom);
 342     if (TMath::Abs(invBenMom)<=1.04) {
 343       AliMUONTrack track;
 344       AliMUONESDInterface::ESDToMUON(*esdTrack, track);
 345       fAlign->ProcessTrack(&track);
 346       fAlign->LocalFit(fTrackOk++,trackParams,0);
 347     }
 348     fTrackTot++;
 349     cout << "Processed " << fTrackTot << " Tracks." << endl;
 350   }
 351
 352   // Post final data. Write histo list to a file with option "RECREATE"
 353   PostData(0,fList);
 354 }
 355
 356 //________________________________________________________________________
 357 void AliMUONAlignmentTask::FinishTaskOutput()
 358 {
 359   /// Called once per task on the client machine at the end of the analysis.
 360
 361   cout << "Processed " << fTrackTot << " Tracks." << endl;
 362   // Perform global fit
 363   fAlign->GlobalFit(fParameters,fErrors,fPulls);
 364
 365   cout << "Done with GlobalFit " << endl;
 366
 367 //   // Update pointers reading them from the output slot
 368 //   fList = (TList*)GetOutputData(0);
 369 //   fMSDEx = (TGraphErrors*)fList->At(0);
 370 //   fMSDEy = (TGraphErrors*)fList->At(1);
 371 //   fMSDEz = (TGraphErrors*)fList->At(2);
 372 //   fMSDEp = (TGraphErrors*)fList->At(3);
 373
 374   // Store results
 375   Double_t deId[156] = {0};
 376   Double_t msdEx[156] = {0};
 377   Double_t msdEy[156] = {0};
 378   Double_t msdEz[156] = {0};
 379   Double_t msdEp[156] = {0};
 380   Double_t deIdErr[156] = {0};
 381   Double_t msdExErr[156] = {0};
 382   Double_t msdEyErr[156] = {0};
 383   Double_t msdEzErr[156] = {0};
 384   Double_t msdEpErr[156] = {0};
 385   Int_t lNDetElem = 4*2+4*2+18*2+26*2+26*2;
 386   Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
 387   // Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
 388   Int_t idOffset = 0; // 400
 389   Int_t lSDetElemCh = 0;
 390   for(Int_t iDE=0; iDE<lNDetElem; iDE++){
 391     deIdErr[iDE] = 0.;
 392     deId[iDE] = idOffset+100;
 393     deId[iDE] += iDE;
 394     lSDetElemCh = 0;
 395     for(Int_t iCh=0; iCh<9; iCh++){
 396       lSDetElemCh += lNDetElemCh[iCh];
 397       if (iDE>=lSDetElemCh) {
 398         deId[iDE] += 100;
 399         deId[iDE] -= lNDetElemCh[iCh];
 400       }
 401     }
 402     msdEx[iDE]=fParameters[3*iDE+0];
 403     msdEy[iDE]=fParameters[3*iDE+1];
 404     msdEz[iDE]=fParameters[3*iDE+3];
 405     msdEp[iDE]=fParameters[3*iDE+2];
 406     msdExErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+0);
 407     msdEyErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+1);
 408     msdEzErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+3);
 409     msdEpErr[iDE]=(Double_t)fAlign->GetParError(3*iDE+2);
 410     fMSDEx->SetPoint(iDE,deId[iDE],fParameters[3*iDE+0]);
 411     fMSDEx->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+0));
 412     fMSDEy->SetPoint(iDE,deId[iDE],fParameters[3*iDE+1]);
 413     fMSDEy->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+1));
 414     fMSDEp->SetPoint(iDE,deId[iDE],fParameters[3*iDE+2]);
 415     fMSDEp->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+2));
 416     fMSDEz->SetPoint(iDE,deId[iDE],fParameters[3*iDE+3]);
 417     fMSDEz->SetPoint(iDE,deIdErr[iDE],(Double_t)fAlign->GetParError(3*iDE+3));
 418   }
 419
 420   // Post final data. Write histo list to a file with option "RECREATE"
 421   PostData(0,fList);
 422
 423   // Re Align
 424   AliMUONGeometryTransformer *newTransform = fAlign->ReAlign(fTransform,fParameters,true);
 425   newTransform->WriteTransformations("transform2ReAlign.dat");
 426
 427   // Generate realigned data in local cdb
 428   const TClonesArray* array = newTransform->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   AliCDBManager* cdbManager = AliCDBManager::Instance();
 435   cdbManager->SetDefaultStorage(fDefaultOCDB.Data());
 436   cdbManager->SetSpecificStorage("MUON/Align/Data",fMisAlignOCDB.Data());
 437
 438   AliCDBMetaData* cdbData = new AliCDBMetaData();
 439   cdbData->SetResponsible("Dimuon Offline project");
 440   cdbData->SetComment("MUON alignment objects with residual misalignment");
 441   AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
 442   cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
 443
 444 }
 445
 446 //________________________________________________________________________
 447 void AliMUONAlignmentTask::Terminate(const Option_t*)
 448 {
 449   /// Called once per task on the client machine at the end of the analysis.
 450
 451 }
 452
 453 //-----------------------------------------------------------------------
 454 void AliMUONAlignmentTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
 455 {
 456   /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
 457
 458   // Load mapping
 459   AliCDBManager* man = AliCDBManager::Instance();
 460   man->SetDefaultStorage(defaultOCDB);
 461   man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
 462   man->Print();
 463   man->SetRun(fLastRunNumber);
 464   if ( ! AliMpCDB::LoadDDLStore() ) {
 465     Error("MUONRefit","Could not access mapping from OCDB !");
 466     exit(-1);
 467   }
 468
 469   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
 470   if (!gGeoManager) {
 471     AliGeomManager::LoadGeometry(geoFilename);
 472     if (!gGeoManager) {
 473       Error("AliMUONReAlignTask", "getting geometry from file %s failed", "generated/galice.root");
 474       return;
 475     }
 476   }
 477
 478   // set mag field
 479   if (!TGeoGlobalMagField::Instance()->GetField()) {
 480     printf("Loading field map...\n");
 481     AliGRPManager *grpMan = new AliGRPManager();
 482     grpMan->ReadGRPEntry();
 483     grpMan->SetMagField();
 484     delete grpMan;
 485   }
 486   // set the magnetic field for track extrapolations
 487   AliMUONTrackExtrap::SetField();
 488
 489 }