MUON/AliMUONReAlignTask.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 AliMUONReAlignTask
  20 /// AliAnalysisTask to realign the MUON spectrometer.
  21 /// The Task reads as input ESDs moves the clusters of a MUONTrack acoording
  22 /// to the re aligned geometry taken from a misalignment file in the OCDB
  23 /// and refit the track. Then it produces a AliMUONClusterInfo object for each
  24 /// cluster. The output is a TTree of AliMUONClusterInfo
  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 <TTree.h>
  34 #include <TChain.h>
  35 #include <TClonesArray.h>
  36 #include <TRandom.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 "AliGeomManager.h"
  50
  51 #include "AliMpConstants.h"
  52 #include "AliMpCDB.h"
  53 #include "AliMpSegmentation.h"
  54 #include "AliMpVSegmentation.h"
  55 #include "AliMpPad.h"
  56 #include "AliMUONCalibrationData.h"
  57 #include "AliMUONVCalibParam.h"
  58 #include "AliMUONPadInfo.h"
  59 #include "AliMUONClusterInfo.h"
  60 #include "AliMUONRecoParam.h"
  61 #include "AliMUONESDInterface.h"
  62 #include "AliMUONRefitter.h"
  63 #include "AliMUONRecoParam.h"
  64 #include "AliMUONVDigit.h"
  65 #include "AliMUONVCluster.h"
  66 #include "AliMUONTrack.h"
  67 #include "AliMUONVTrackStore.h"
  68 #include "AliMUONVDigitStore.h"
  69 #include "AliMUONTrackParam.h"
  70 #include "AliMUONTrackExtrap.h"
  71 #include "AliMUONGeometryTransformer.h"
  72
  73 #include "AliMUONReAlignTask.h"
  74
  75 ///\cond CLASSIMP
  76 ClassImp(AliMUONReAlignTask)
  77 ///\endcond
  78
  79 //________________________________________________________________________
  80 AliMUONReAlignTask::AliMUONReAlignTask(const char *name, const char *geofilename, const char *defaultocdb, const char *misalignocdb)
  81   : AliAnalysisTask(name, ""),
  82     fESD(0x0),
  83     fClusterInfoTree(0x0),
  84     fClusterInfo(0x0),
  85     fESDInterface(0x0),
  86     fRefitter(0x0),
  87     fRecoParam(0x0),
  88     fGeoFilename(geofilename),
  89     fMisAlignOCDB(misalignocdb),
  90     fDefaultOCDB(defaultocdb),
  91     fGeoTransformer(0x0),
  92     fNewGeoTransformer(0x0),
  93     fGainStore(0x0),
  94     fPedStore(0x0),
  95     fPrintLevel(0),
  96     fLastRun(-1)
  97 {
  98   /// Default Constructor
  99   // Define input and output slots here
 100   // Input slot #0 works with a TChain
 101   DefineInput(0, TChain::Class());
 102   // Output slot #0 writes a TTree
 103   DefineOutput(0, TTree::Class());
 104
 105   fClusterInfo = new AliMUONClusterInfo();
 106   fESDInterface = new AliMUONESDInterface();
 107   fGeoTransformer = new AliMUONGeometryTransformer();
 108   fNewGeoTransformer = new AliMUONGeometryTransformer();
 109 }
 110
 111 //________________________________________________________________________
 112 AliMUONReAlignTask::AliMUONReAlignTask(const AliMUONReAlignTask& obj)
 113   : AliAnalysisTask(obj),
 114     fESD(0x0),
 115     fClusterInfoTree(0x0),
 116     fClusterInfo(0x0),
 117     fESDInterface(0x0),
 118     fRefitter(0x0),
 119     fRecoParam(0x0),
 120     fGeoFilename(""),
 121     fMisAlignOCDB(""),
 122     fDefaultOCDB(""),
 123     fGeoTransformer(0x0),
 124     fNewGeoTransformer(0x0),
 125     fGainStore(0x0),
 126     fPedStore(0x0),
 127     fPrintLevel(0),
 128     fLastRun(-1)
 129 {
 130   /// Copy constructor
 131   fESD = obj.fESD;
 132   fClusterInfoTree = obj.fClusterInfoTree;
 133   fClusterInfo = obj.fClusterInfo;
 134   fESDInterface = obj.fESDInterface;
 135   fRefitter = obj.fRefitter;
 136   fRecoParam = obj.fRecoParam;
 137   fGeoFilename = obj.fGeoFilename;
 138   fMisAlignOCDB = obj.fMisAlignOCDB;
 139   fDefaultOCDB = obj.fDefaultOCDB;
 140   fGeoTransformer = obj.fGeoTransformer;
 141   fNewGeoTransformer = obj.fNewGeoTransformer;
 142   fGainStore = obj.fGainStore;
 143   fPedStore = obj.fPedStore;
 144   fPrintLevel = obj.fPrintLevel;
 145   fLastRun = obj.fLastRun;
 146 }
 147
 148 //________________________________________________________________________________
 149 AliMUONReAlignTask& AliMUONReAlignTask::operator=(const AliMUONReAlignTask& other)
 150 {
 151   /// Assignment
 152   if(&other == this) return *this;
 153   AliAnalysisTask::operator=(other);
 154   fESD = other.fESD;
 155   fClusterInfoTree = other.fClusterInfoTree;
 156   fClusterInfo = other.fClusterInfo;
 157   fESDInterface = other.fESDInterface;
 158   fRefitter = other.fRefitter;
 159   fRecoParam = other.fRecoParam;
 160   fGeoFilename = other.fGeoFilename;
 161   fMisAlignOCDB = other.fMisAlignOCDB;
 162   fDefaultOCDB = other.fDefaultOCDB;
 163   fGeoTransformer = other.fGeoTransformer;
 164   fNewGeoTransformer = other.fNewGeoTransformer;
 165   fGainStore = other.fGainStore;
 166   fPedStore = other.fPedStore;
 167   fPrintLevel = other.fPrintLevel;
 168   fLastRun = other.fLastRun;
 169
 170   return *this;
 171 }
 172
 173
 174 //________________________________________________________________________
 175 AliMUONReAlignTask::~AliMUONReAlignTask()
 176 {
 177   /// Destructor
 178   if (fESDInterface) delete fESDInterface;
 179   if (fGeoTransformer) delete fGeoTransformer;
 180   if (fNewGeoTransformer) delete fNewGeoTransformer;
 181 }
 182
 183 //________________________________________________________________________
 184 void AliMUONReAlignTask::LocalInit()
 185 {
 186   /// Local initialization, called once per task on the client machine
 187   /// where the analysis train is assembled
 188   AliMpCDB::LoadMpSegmentation();
 189
 190   // prepare the refitting
 191   gRandom->SetSeed(0);
 192   Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
 193
 194   fRefitter = new AliMUONRefitter(fRecoParam);
 195   fRefitter->Connect(fESDInterface);
 196
 197   // Original geotransformer
 198   fGeoTransformer->LoadGeometryData();
 199   // Apply mis alignments
 200   AliGeomManager::ApplyAlignObjsFromCDB("MUON");
 201   // New geotransformer
 202   fNewGeoTransformer->LoadGeometryData();
 203
 204 }
 205
 206 //________________________________________________________________________
 207 void AliMUONReAlignTask::ConnectInputData(Option_t *)
 208 {
 209   /// Connect ESD here. Called on each input data change.
 210   // Connect ESD here
 211   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
 212   if (!esdTree) {
 213     Printf("ERROR: Could not read chain from input slot 0");
 214   }
 215   else {
 216     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 217     if (!esdH) {
 218       Printf("ERROR: Could not get ESDInputHandler");
 219     }
 220     else {
 221       fESD = esdH->GetEvent();
 222     }
 223   }
 224 }
 225
 226 //________________________________________________________________________
 227 void AliMUONReAlignTask::CreateOutputObjects()
 228 {
 229   /// Executed once on each worker (machine actually running the analysis code)
 230   //
 231   // This method has to be called INSIDE the user redefined CreateOutputObjects
 232   // method, before creating each object corresponding to the output containers
 233   // that are to be written to a file. This need to be done in general for the big output
 234   // objects that may not fit memory during processing.
 235   OpenFile(0);
 236
 237   fClusterInfoTree = new TTree("clusterInfoTree","clusterInfoTree");
 238   fClusterInfoTree->Branch("clusterInfo", &fClusterInfo, 32000, 99);
 239
 240 }
 241
 242 //________________________________________________________________________
 243 void AliMUONReAlignTask::Exec(Option_t *)
 244 {
 245   /// Main loop, called for each event
 246   if (!fESD) {
 247     Printf("ERROR: fESD not available");
 248     return;
 249   }
 250
 251   // Prepare Gain and Pedestal store
 252   if (!fGainStore || fLastRun!=fESD->GetRunNumber()){
 253     fGainStore = AliMUONCalibrationData::CreateGains(fESD->GetRunNumber());
 254     fLastRun = fESD->GetRunNumber();
 255   }
 256   if (!fPedStore || fLastRun!=fESD->GetRunNumber()){
 257     fPedStore = AliMUONCalibrationData::CreatePedestals(fESD->GetRunNumber());
 258     fLastRun = fESD->GetRunNumber();
 259   }
 260
 261   AliMUONPadInfo padInfo;
 262
 263   Int_t nTracks = (Int_t)fESD->GetNumberOfMuonTracks();
 264   if (nTracks < 1) return;
 265
 266   // load the current event
 267   fESDInterface->LoadEvent(*fESD);
 268   AliMUONVDigitStore* digitStore = fESDInterface->GetDigits();
 269
 270   Double_t lX = 0.;
 271   Double_t lY = 0.;
 272   Double_t lZ = 0.;
 273   Double_t gX = 0.;
 274   Double_t gY = 0.;
 275   Double_t gZ = 0.;
 276   // loop over cluster to modify their position
 277   AliMUONVCluster *cluster;
 278   TIter next(fESDInterface->CreateClusterIterator());
 279   while ((cluster = static_cast<AliMUONVCluster*>(next()))) {
 280     //       cout << "Original cluster" << endl;
 281     //       cluster->Print();
 282     gX = cluster->GetX();
 283     gY = cluster->GetY();
 284     gZ = cluster->GetZ();
 285     fGeoTransformer->Global2Local(cluster->GetDetElemId(),gX,gY,gZ,lX,lY,lZ);
 286     fNewGeoTransformer->Local2Global(cluster->GetDetElemId(),lX,lY,lZ,gX,gY,gZ);
 287     cluster->SetXYZ(gX,gY,gZ);
 288     //       cout << "Aligned cluster" << endl;
 289     //       cluster->Print();
 290   }
 291
 292   // refit the tracks from digits
 293   AliMUONVTrackStore* newTrackStore = fRefitter->ReconstructFromClusters();
 294
 295   //----------------------------------------------//
 296   // ------ fill new ESD and print results ------ //
 297   //----------------------------------------------//
 298   // loop over the list of ESD tracks
 299   TClonesArray *esdTracks = (TClonesArray*) fESD->FindListObject("MuonTracks");
 300   for (Int_t iTrack = 0; iTrack <  nTracks; iTrack++) {
 301     // get the ESD track
 302     AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
 303
 304     // skip ghost tracks (leave them unchanged in the new ESD file)
 305     if (!esdTrack->ContainTrackerData()) continue;
 306
 307     // get the corresponding MUON track
 308     AliMUONTrack* track = fESDInterface->FindTrack(esdTrack->GetUniqueID());
 309
 310     // Find the corresponding re-fitted MUON track
 311     AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
 312
 313 //     // replace the content of the current ESD track or remove it if the refitting has failed
 314 //     if (newTrack) {
 315 //       Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
 316 //       AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex, esdInterface.GetDigits());
 317 //       } else {
 318 //       esdTracks->Remove(esdTrack);
 319 //       }
 320
 321     // print initial and re-fitted track parameters at first cluster if any
 322     if (fPrintLevel>0) {
 323       cout<<"            ----------------track #"<<iTrack+1<<"----------------"<<endl;
 324       cout<<"before refit:"<<endl;
 325       AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
 326       param->Print("FULL");
 327       if (fPrintLevel>1) param->GetCovariances().Print();
 328       if (!newTrack) continue;
 329       cout<<"after refit:"<<endl;
 330       param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
 331       param->Print("FULL");
 332       if (fPrintLevel>1) param->GetCovariances().Print();
 333       cout<<"            ----------------------------------------"<<endl;
 334     }
 335
 336     // Cluster Info part
 337     UInt_t muonClusterMap = BuildClusterMap(*newTrack);
 338
 339     // loop over clusters
 340     AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->First());
 341     while (trackParam) {
 342       fClusterInfo->Clear("C");
 343
 344       // fill cluster info
 345       cluster = trackParam->GetClusterPtr();
 346       fClusterInfo->SetRunId(fESD->GetRunNumber());
 347       fClusterInfo->SetEventId(fESD->GetEventNumberInFile());
 348       fClusterInfo->SetZ(cluster->GetZ());
 349       fClusterInfo->SetClusterId(cluster->GetUniqueID());
 350       fClusterInfo->SetClusterXY(cluster->GetX(), cluster->GetY());
 351       fClusterInfo->SetClusterXYErr(cluster->GetErrX(), cluster->GetErrY());
 352       fClusterInfo->SetClusterChi2(cluster->GetChi2());
 353       fClusterInfo->SetClusterCharge(cluster->GetCharge());
 354
 355       // fill track info
 356       fClusterInfo->SetTrackId(newTrack->GetUniqueID());
 357       fClusterInfo->SetTrackXY(trackParam->GetNonBendingCoor(), trackParam->GetBendingCoor());
 358       fClusterInfo->SetTrackThetaXY(TMath::ATan(trackParam->GetNonBendingSlope()), TMath::ATan(trackParam->GetBendingSlope()));
 359       fClusterInfo->SetTrackP(trackParam->P());
 360       const TMatrixD paramCov = trackParam->GetCovariances();
 361       fClusterInfo->SetTrackXYErr(TMath::Sqrt(paramCov(0,0)), TMath::Sqrt(paramCov(2,2)));
 362       fClusterInfo->SetTrackChi2(newTrack->GetNormalizedChi2());
 363       fClusterInfo->SetTrackCharge((Short_t)trackParam->GetCharge());
 364       fClusterInfo->SetTrackNHits(newTrack->GetNClusters());
 365       fClusterInfo->SetTrackChamberHitMap(muonClusterMap);
 366
 367       // fill pad info if available
 368       for (Int_t i=0; i<cluster->GetNDigits(); i++) {
 369         AliMUONVDigit* digit = digitStore->FindObject(cluster->GetDigitId(i));
 370         if (!digit) continue;
 371
 372         // pad location
 373         const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
 374           GetMpSegmentation(digit->DetElemId(),AliMp::GetCathodType(digit->Cathode()));
 375         AliMpPad pad = seg->PadByIndices(digit->PadX(), digit->PadY());
 376
 377         // calibration parameters
 378         AliMUONVCalibParam* ped =  fPedStore ? static_cast<AliMUONVCalibParam*>(fPedStore->FindObject(digit->DetElemId(), digit->ManuId())) : 0x0;
 379         AliMUONVCalibParam* gain = fGainStore ? static_cast<AliMUONVCalibParam*>(fGainStore->FindObject(digit->DetElemId(), digit->ManuId())) : 0x0;
 380         Int_t manuChannel = digit->ManuChannel();
 381         Int_t planeType = 0;
 382         if ( digit->ManuId() & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) {
 383           planeType = 1;
 384         }
 385
 386         // fill pad info
 387         padInfo.SetPadId(digit->GetUniqueID());
 388         padInfo.SetPadPlaneType(planeType);
 389         padInfo.SetPadXY(pad.GetPositionX(), pad.GetPositionY());
 390         padInfo.SetPadDimXY(pad.GetDimensionX(), pad.GetDimensionY());
 391         padInfo.SetPadCharge((Double_t)digit->Charge());
 392         padInfo.SetPadADC(digit->ADC());
 393         padInfo.SetSaturated(digit->IsSaturated());
 394         padInfo.SetCalibrated(digit->IsCalibrated());
 395         if (ped) {
 396           padInfo.SetPedestal(ped->ValueAsFloatFast(manuChannel,0), ped->ValueAsFloatFast(manuChannel,1));
 397         } else {
 398           padInfo.SetPedestal(-250.,-5.);
 399         }
 400         if (gain) {
 401           padInfo.SetGain(gain->ValueAsFloatFast(manuChannel,0), gain->ValueAsFloatFast(manuChannel,1),
 402                           gain->ValueAsIntFast(manuChannel,2), gain->ValueAsIntFast(manuChannel,3));
 403         } else {
 404           padInfo.SetGain(-1.,-0.1,-4095,-1);
 405         }
 406         fClusterInfo->AddPad(padInfo);
 407       }
 408
 409       // fill cluster info tree
 410       fClusterInfoTree->Fill();
 411       trackParam = static_cast<AliMUONTrackParam*>(track->GetTrackParamAtCluster()->After(trackParam));
 412     }
 413   }
 414   // free memory
 415   delete newTrackStore;
 416
 417
 418   // Post final data. Write histo list to a file with option "RECREATE"
 419   PostData(0,fClusterInfoTree);
 420 }
 421
 422 //________________________________________________________________________
 423 void AliMUONReAlignTask::Terminate(const Option_t*)
 424 {
 425   /// Called once per task on the client machine at the end of the analysis.
 426
 427 }
 428
 429 //-----------------------------------------------------------------------
 430 void AliMUONReAlignTask::Prepare(const char* geoFilename, const char* defaultOCDB, const char* misAlignOCDB)
 431 {
 432   /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
 433
 434   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
 435   if (!gGeoManager) {
 436     AliGeomManager::LoadGeometry(geoFilename);
 437     if (!gGeoManager) {
 438       Error("AliMUONReAlignTask", "getting geometry from file %s failed", "generated/galice.root");
 439       return;
 440     }
 441   }
 442
 443   // Load mapping
 444   AliCDBManager* man = AliCDBManager::Instance();
 445   man->SetDefaultStorage(defaultOCDB);
 446   man->SetSpecificStorage("MUON/Align/Data",misAlignOCDB);
 447   man->Print();
 448   man->SetRun(0);
 449   if ( ! AliMpCDB::LoadDDLStore() ) {
 450     Error("MUONRefit","Could not access mapping from OCDB !");
 451     exit(-1);
 452   }
 453
 454   // set mag field
 455   if (!TGeoGlobalMagField::Instance()->GetField()) {
 456     printf("Loading field map...\n");
 457     AliGRPManager *grpMan = new AliGRPManager();
 458     grpMan->ReadGRPEntry();
 459     grpMan->SetMagField();
 460     delete grpMan;
 461   }
 462   // set the magnetic field for track extrapolations
 463   AliMUONTrackExtrap::SetField();
 464
 465   // Load initial reconstruction parameters from OCDB
 466   // reconstruction parameters
 467   fRecoParam = AliMUONRecoParam::GetCosmicParam();
 468
 469   // digit selection
 470   fRecoParam->SetPadGoodnessMask(0x400BE80);
 471 //   TString caliboption = caliboption1;
 472 //   if ( calib == 2 ) caliboption = caliboption2;
 473 //   fRecoParam->SetCalibrationMode("NOGAIN"caliboption.Data());
 474   fRecoParam->SetCalibrationMode("NOGAIN");
 475
 476   // chamber resolution (incuding misalignment)
 477   for (Int_t iCh=0; iCh<10; iCh++) {
 478     fRecoParam->SetDefaultNonBendingReso(iCh,0.4);
 479     fRecoParam->SetDefaultBendingReso(iCh,0.4);
 480   }
 481   fRecoParam->SetMaxNonBendingDistanceToTrack(10.);
 482   fRecoParam->SetMaxBendingDistanceToTrack(10.);
 483
 484   // cut on (non)bending slopes
 485   //fRecoParam->SetMaxNonBendingSlope(0.6);
 486   //fRecoParam->SetMaxBendingSlope(0.6);
 487
 488   // tracking algorithm
 489   //  fRecoParam->MakeMoreTrackCandidates(kTRUE);
 490   fRecoParam->RequestStation(0, kFALSE);
 491   fRecoParam->RequestStation(2, kFALSE);
 492 //   fRecoParam->RequestStation(3, kFALSE);
 493 //   fRecoParam->RequestStation(4, kFALSE);
 494   fRecoParam->SetSigmaCutForTracking(7.);
 495   fRecoParam->ImproveTracks(kTRUE, 7.);
 496   Info("MUONRefit", "\n initial recontruction parameters:");
 497   fRecoParam->Print("FULL");
 498
 499   AliMUONESDInterface::ResetTracker(fRecoParam);
 500 }
 501
 502 //-----------------------------------------------------------------------
 503 UInt_t AliMUONReAlignTask::BuildClusterMap(AliMUONTrack &track)
 504 {
 505   /// Build the map of clusters in tracking chambers
 506
 507   UInt_t muonClusterMap = 0;
 508
 509   AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
 510   while (trackParam) {
 511
 512     muonClusterMap |= BIT(trackParam->GetClusterPtr()->GetChamberId());
 513
 514     trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
 515   }
 516
 517   return muonClusterMap;
 518
 519 }