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