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