MUON/AliMUONReAlignTask.cxx

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