MUON/AliMUONAlignmentTask.cxx

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