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