PWG/muondep/AliAnalysisTaskMuonRefit.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 // ROOT includes
  19 #include <TString.h>
  20 #include <TList.h>
  21 #include <TGeoManager.h>
  22 #include <TGeoGlobalMagField.h>
  23 #include <TPRegexp.h>
  24
  25 // STEER includes
  26 #include "AliESDEvent.h"
  27 #include "AliESDMuonTrack.h"
  28 #include "AliCDBManager.h"
  29 #include "AliGRPManager.h"
  30 #include "AliGRPObject.h"
  31 #include "AliCDBStorage.h"
  32 #include "AliGeomManager.h"
  33 #include "AliMagF.h"
  34
  35 // ANALYSIS includes
  36 #include "AliAnalysisManager.h"
  37 #include "AliInputEventHandler.h"
  38 #include "AliAnalysisTaskMuonRefit.h"
  39
  40 // MUON includes
  41 #include "AliMUONCDB.h"
  42 #include "AliMUONConstants.h"
  43 #include "AliMUONRecoParam.h"
  44 #include "AliMUONESDInterface.h"
  45 #include "AliMUONRefitter.h"
  46 #include "AliMUONTrack.h"
  47 #include "AliMUONTrackParam.h"
  48 #include "AliMUONTrackExtrap.h"
  49 #include "AliMUONVCluster.h"
  50 #include "AliMUONVDigit.h"
  51 #include "AliMUONVDigitStore.h"
  52 #include "AliMUONVTrackStore.h"
  53 #include "AliMUONLocalTrigger.h"
  54 #include "AliMUONTriggerTrack.h"
  55 #include "AliMUONVTriggerStore.h"
  56 #include "AliMUONVTriggerTrackStore.h"
  57 #include "AliMUONGeometryTransformer.h"
  58 #include "AliMUONTriggerCircuit.h"
  59 #include "AliMUONTrackHitPattern.h"
  60 #include "AliMUONVTrackReconstructor.h"
  61
  62 // MUON mapping includes
  63 #include "AliMpSegmentation.h"
  64 #include "AliMpVSegmentation.h"
  65 #include "AliMpConstants.h"
  66 #include "AliMpDDLStore.h"
  67 #include "AliMpManuStore.h"
  68 #include "AliMpPad.h"
  69 #include "AliMpDetElement.h"
  70 #include "AliMpCathodType.h"
  71
  72 #ifndef SafeDelete
  73 #define SafeDelete(x) if (x != NULL) { delete x; x = NULL; }
  74 #endif
  75
  76 ClassImp(AliAnalysisTaskMuonRefit)
  77
  78 //________________________________________________________________________
  79 AliAnalysisTaskMuonRefit::AliAnalysisTaskMuonRefit() :
  80 AliAnalysisTaskSE(),
  81 fDefaultStorage(""),
  82 fImproveTracks(kFALSE),
  83 fSigmaCut(-1.),
  84 fSigmaCutForTrigger(-1.),
  85 fReAlign(kFALSE),
  86 fOldAlignStorage(""),
  87 fNewAlignStorage(""),
  88 fOldGeoTransformer(NULL),
  89 fNewGeoTransformer(NULL),
  90 fField(""),
  91 fRemoveMonoCathCl(kFALSE),
  92 fCheckAllPads(kFALSE),
  93 fTagBadTracks(kFALSE),
  94 fKeepOldParam(kFALSE),
  95 fTriggerCircuit(NULL),
  96 fESDInterface(NULL),
  97 fRefitter(NULL),
  98 fTrackHitPattern(NULL)
  99 {
 100   /// Default constructor
 101   for (Int_t i = 0; i < 10; i++) ResetClusterResolution(i, -1., -1.);
 102 }
 103
 104 //________________________________________________________________________
 105 AliAnalysisTaskMuonRefit::AliAnalysisTaskMuonRefit(const char *name) :
 106 AliAnalysisTaskSE(name),
 107 fDefaultStorage("raw://"),
 108 fImproveTracks(kFALSE),
 109 fSigmaCut(-1.),
 110 fSigmaCutForTrigger(-1.),
 111 fReAlign(kFALSE),
 112 fOldAlignStorage(""),
 113 fNewAlignStorage(""),
 114 fOldGeoTransformer(NULL),
 115 fNewGeoTransformer(NULL),
 116 fField(""),
 117 fRemoveMonoCathCl(kFALSE),
 118 fCheckAllPads(kFALSE),
 119 fTagBadTracks(kFALSE),
 120 fKeepOldParam(kFALSE),
 121 fTriggerCircuit(NULL),
 122 fESDInterface(NULL),
 123 fRefitter(NULL),
 124 fTrackHitPattern(NULL)
 125 {
 126   /// Constructor
 127   fBranchNames = "ESD:AliESDRun.,AliESDHeader.,MuonTracks,MuonClusters,MuonPads";
 128   for (Int_t i = 0; i < 10; i++) ResetClusterResolution(i, -1., -1.);
 129 }
 130
 131 //________________________________________________________________________
 132 AliAnalysisTaskMuonRefit::~AliAnalysisTaskMuonRefit()
 133 {
 134   /// Destructor
 135   SafeDelete(fOldGeoTransformer);
 136   SafeDelete(fNewGeoTransformer);
 137   SafeDelete(fTriggerCircuit);
 138   SafeDelete(fESDInterface);
 139   SafeDelete(fRefitter);
 140   SafeDelete(fTrackHitPattern);
 141 }
 142
 143 //___________________________________________________________________________
 144 void AliAnalysisTaskMuonRefit::UserCreateOutputObjects()
 145 {
 146 }
 147
 148 //________________________________________________________________________
 149 void AliAnalysisTaskMuonRefit::UserExec(Option_t *)
 150 {
 151   /// Main event loop
 152
 153   // check if refitter properly created (i.e. OCDB data properly loaded)
 154   if (!fRefitter) AliFatal("Problem occur while loading OCDB objects");
 155
 156   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 157   if (!esd) return;
 158
 159   LoadBranches();
 160
 161   Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
 162   if (nTracks < 1) return;
 163
 164   TList oldGhosts;
 165   oldGhosts.SetOwner(kFALSE);
 166   UInt_t firstGhostId = 0xFFFFFFFF - 1;
 167
 168   TList esdTracksToRemove;
 169   esdTracksToRemove.SetOwner(kFALSE);
 170
 171   // load the current event
 172   fESDInterface->LoadEvent(*esd, kFALSE);
 173
 174   // remove clusters from ESD (keep digits as they will not change, just eventually not used anymore)
 175   esd->FindListObject("MuonClusters")->Clear("C");
 176
 177   // modify clusters
 178   Int_t firstClusterIndex = 0;
 179   AliMUONVCluster* cluster = 0x0;
 180   TIter nextCluster(fESDInterface->CreateClusterIterator());
 181   while ((cluster = static_cast<AliMUONVCluster*>(nextCluster()))) {
 182     Int_t clIndex = cluster->GetClusterIndex(cluster->GetUniqueID());
 183     if (clIndex >= firstClusterIndex) firstClusterIndex = clIndex+1;
 184     ModifyCluster(*cluster);
 185   }
 186
 187   // to do not mix with old clusters in case we re-clusterize
 188   fRefitter->SetFirstClusterIndex(firstClusterIndex);
 189
 190   // refit the tracks from clusters
 191   AliMUONVTrackStore* newTrackStore = fRefitter->ReconstructFromClusters();
 192
 193   // reset the trigger part
 194   TIter nextNewTrack(newTrackStore->CreateIterator());
 195   AliMUONTrack *newTrack = 0x0;
 196   while ((newTrack = static_cast<AliMUONTrack*>(nextNewTrack()))) {
 197     newTrack->SetMatchTrigger(0);
 198     newTrack->SetLocalTrigger(0,0,0,0,0,0,0);
 199     newTrack->SetChi2MatchTrigger(0.);
 200     newTrack->SetHitsPatternInTrigCh(0);
 201     newTrack->SetHitsPatternInTrigChTrk(0);
 202   }
 203
 204   // reconstruct trigger tracks
 205   AliMUONVTriggerStore *trigStore = fESDInterface->GetTriggers();
 206   AliMUONVTriggerTrackStore *trigTrackStore = AliMUONESDInterface::NewTriggerTrackStore();
 207   if (!trigTrackStore) return;
 208   AliMUONVTrackReconstructor* tracker = AliMUONESDInterface::GetTracker();
 209   tracker->EventReconstructTrigger(*fTriggerCircuit, *trigStore, *trigTrackStore);
 210
 211   // recover the hit pattern for all trigger tracks
 212   TClonesArray *esdTracks = (TClonesArray*) esd->FindListObject("MuonTracks");
 213   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 214     AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
 215     if (!esdTrack->ContainTriggerData()) continue;
 216     // use the UniqueID of the local trigger in case several tracks match the same trigger
 217     AliMUONLocalTrigger *locTrg = fESDInterface->FindLocalTrigger(esdTrack->LoCircuit());
 218     AliMUONTriggerTrack *trigTrack = (AliMUONTriggerTrack*) trigTrackStore->FindObject(locTrg->GetUniqueID());
 219     trigTrack->SetHitsPatternInTrigCh(esdTrack->GetHitsPatternInTrigCh());
 220   }
 221
 222   // match tracker/trigger tracks
 223   const Int_t kFirstTrigCh = AliMUONConstants::NTrackingCh();
 224   const Float_t kZFilterOut = AliMUONConstants::MuonFilterZEnd();
 225   const Float_t kFilterThickness = kZFilterOut-AliMUONConstants::MuonFilterZBeg();
 226   nextNewTrack.Reset();
 227   while ((newTrack = static_cast<AliMUONTrack*>(nextNewTrack()))) {
 228     AliMUONTrackParam trackParam(*((AliMUONTrackParam*) (newTrack->GetTrackParamAtCluster()->Last())));
 229     AliMUONTrackExtrap::ExtrapToZCov(&trackParam, kZFilterOut);
 230     AliMUONTrackExtrap::AddMCSEffect(&trackParam, kFilterThickness, AliMUONConstants::MuonFilterX0());
 231     AliMUONTrackExtrap::ExtrapToZCov(&trackParam, AliMUONConstants::DefaultChamberZ(kFirstTrigCh));
 232     AliMUONTriggerTrack *matchedTriggerTrack = fTrackHitPattern->MatchTriggerTrack(newTrack, trackParam, *trigTrackStore, *trigStore);
 233     if ( matchedTriggerTrack ) newTrack->SetHitsPatternInTrigCh(matchedTriggerTrack->GetHitsPatternInTrigCh());
 234   }
 235
 236   // loop over the list of ESD tracks
 237   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 238
 239     // get the ESD track
 240     AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
 241
 242     // keep the memory of old ghosts to concerve their Id if their are still ghosts after refitting
 243     // and remember to remove them before adding the new ghosts
 244     if (!esdTrack->ContainTrackerData()) {
 245       AliMUONLocalTrigger *locTrg = (esdTrack->ContainTriggerData()) ? fESDInterface->FindLocalTrigger(esdTrack->LoCircuit()) : 0x0;
 246       if (locTrg && !locTrg->IsNull()) oldGhosts.AddLast(locTrg);
 247       if (esdTrack->GetUniqueID() <= firstGhostId) firstGhostId = esdTrack->GetUniqueID()-1;
 248       esdTracksToRemove.AddLast(esdTrack);
 249       continue;
 250     }
 251
 252     // Find the corresponding re-fitted MUON track
 253     newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
 254
 255     // replace the content of the current ESD track or remove it
 256     if (newTrack && (!fImproveTracks || newTrack->IsImproved())) {
 257
 258       // eventually tag the tracks which do not match the trigger while they were before
 259       Bool_t noLongerMatched = (fTagBadTracks && esdTrack->ContainTriggerData() && newTrack->GetMatchTrigger() <= 0);
 260
 261       // find the corresponding trigger info and eventually restore the old one
 262       AliMUONLocalTrigger *locTrg = 0x0;
 263       if (newTrack->GetMatchTrigger() > 0) locTrg = fESDInterface->FindLocalTrigger(newTrack->LoCircuit());
 264       else if (noLongerMatched && !fKeepOldParam) {
 265         locTrg = fESDInterface->FindLocalTrigger(esdTrack->LoCircuit());
 266         newTrack->SetMatchTrigger(esdTrack->GetMatchTrigger());
 267         newTrack->SetChi2MatchTrigger(esdTrack->GetChi2MatchTrigger());
 268         newTrack->SetHitsPatternInTrigCh(esdTrack->GetHitsPatternInTrigCh());
 269         newTrack->SetHitsPatternInTrigChTrk(esdTrack->GetHitsPatternInTrigChTrk());
 270         newTrack->SetLocalTrigger(esdTrack->LoCircuit(), esdTrack->LoStripX(), esdTrack->LoStripY(),
 271                                   esdTrack->LoDev(), esdTrack->LoLpt(), esdTrack->LoHpt(),
 272                                   esdTrack->GetTriggerWithoutChamber());
 273       }
 274       if (locTrg && locTrg->IsNull()) locTrg = 0x0;
 275
 276       // fill the new track/cluster info
 277       if (!(noLongerMatched && fKeepOldParam)) {
 278
 279         // fill the track info
 280         Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
 281         AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex, locTrg);
 282
 283         // add the clusters if not already there
 284         for (Int_t i = 0; i < newTrack->GetNClusters(); i++) {
 285           AliMUONVCluster *cl = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->UncheckedAt(i))->GetClusterPtr();
 286           if (esd->FindMuonCluster(cl->GetUniqueID())) continue;
 287           AliESDMuonCluster *esdCl = esd->NewMuonCluster();
 288           AliMUONESDInterface::MUONToESD(*cl, *esdCl, kTRUE);
 289         }
 290
 291       } else {
 292
 293         // restore the old clusters if not already done
 294         for (Int_t i = 0; i < esdTrack->GetNClusters(); i++) {
 295           UInt_t clId = esdTrack->GetClusterId(i);
 296           if (esd->FindMuonCluster(clId)) continue;
 297           AliMUONVCluster *cl = fESDInterface->FindCluster(esdTrack->GetUniqueID(), clId);
 298           if (!cl) continue;
 299           AliESDMuonCluster *esdCl = esd->NewMuonCluster();
 300           AliMUONESDInterface::MUONToESD(*cl, *esdCl, kTRUE);
 301         }
 302
 303       }
 304
 305       // eventually tag the trigger part as bad
 306       esdTrack->ResetBit(BIT(20));
 307       if (noLongerMatched) esdTrack->SetBit(BIT(21));
 308       else esdTrack->ResetBit(BIT(21));
 309
 310     } else {
 311
 312       // simply tag it or remember to remove that track (cannot remove it now as it will create a hole which
 313       // will eventually produce a crash in parallel track loop (AliESDEvent::MoveMuonObjects()))
 314       if (fTagBadTracks) {
 315         esdTrack->SetBit(BIT(20));
 316         if (esdTrack->ContainTriggerData()) esdTrack->SetBit(BIT(21));
 317         else esdTrack->ResetBit(BIT(21));
 318       }
 319       else esdTracksToRemove.AddLast(esdTrack);
 320
 321     }
 322
 323   }
 324
 325   // free memory
 326   delete newTrackStore;
 327
 328   // remove tracks to remove and compress the array of ESD tracks
 329   TIter nextTrackToRemove(&esdTracksToRemove);
 330   AliESDMuonTrack* esdTrack = 0x0;
 331   while ((esdTrack = static_cast<AliESDMuonTrack*>(nextTrackToRemove()))) esdTracks->Remove(esdTrack);
 332   esdTracks->Compress();
 333
 334   // add new ghosts (ignoring tracks marked as bad or no longer matched)
 335   TIter nextGhost(trigStore->CreateIterator());
 336   AliMUONLocalTrigger *locTrg = 0x0;
 337   while ((locTrg = static_cast<AliMUONLocalTrigger*>(nextGhost()))) {
 338     Bool_t alreadyThere = kFALSE;
 339     for (Int_t iTrack = 0; iTrack < esdTracks->GetEntriesFast(); iTrack++) {
 340       esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
 341       alreadyThere = (esdTrack->LoCircuit() == locTrg->LoCircuit() && !esdTrack->TestBit(BIT(21)));
 342       if (alreadyThere) break;
 343     }
 344     if (!alreadyThere) {
 345       UInt_t ghostId = oldGhosts.Contains(locTrg) ? locTrg->GetUniqueID() : firstGhostId--;
 346       AliMUONTriggerTrack *trigTrack = (AliMUONTriggerTrack*) trigTrackStore->FindObject(locTrg->GetUniqueID());
 347       AliMUONESDInterface::MUONToESD(*locTrg, *esd, ghostId, trigTrack);
 348     }
 349   }
 350
 351 }
 352
 353 //________________________________________________________________________
 354 void AliAnalysisTaskMuonRefit::NotifyRun()
 355 {
 356   /// load necessary data from OCDB and create the refitter
 357
 358   // do it only once
 359   if (fRefitter) return;
 360
 361   // set OCDB location
 362   AliCDBManager* cdbm = AliCDBManager::Instance();
 363   if (cdbm->IsDefaultStorageSet()) printf("MuonRefit: CDB default storage already set!\n");
 364   else {
 365     cdbm->SetDefaultStorage(fDefaultStorage.Data());
 366     if (fOldAlignStorage != "none" && !fOldAlignStorage.IsNull())
 367       cdbm->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
 368   }
 369   if (cdbm->GetRun() > -1) printf("MuonRefit: run number already set!\n");
 370   else cdbm->SetRun(fCurrentRunNumber);
 371
 372   // load magnetic field or create it for track extrapolation
 373   if (!TGeoGlobalMagField::Instance()->GetField()) {
 374     if (!fField.IsNull()) {
 375       if (!SetMagField()) return;
 376     } else if (!AliMUONCDB::LoadField()) return;
 377   }
 378
 379   // load mapping
 380   if (!AliMpDDLStore::Instance(kFALSE) || !AliMpManuStore::Instance(kFALSE)) {
 381     if (!AliMUONCDB::LoadMapping()) return;
 382   }
 383
 384   // load recoParam for refitting
 385   AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
 386   if (!recoParam) return;
 387   if (!AliMUONESDInterface::GetTracker()) AliMUONESDInterface::ResetTracker(recoParam);
 388
 389   if (fImproveTracks) {
 390     if (fSigmaCut > 0.) recoParam->ImproveTracks(kTRUE, fSigmaCut);
 391     else recoParam->ImproveTracks(kTRUE);
 392   } else recoParam->ImproveTracks(kFALSE);
 393
 394   if (fSigmaCutForTrigger > 0.) recoParam->SetSigmaCutForTrigger(fSigmaCutForTrigger);
 395   else fSigmaCutForTrigger = recoParam->GetSigmaCutForTrigger();
 396
 397   for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) {
 398     if (fClusterResNB[i] < 0.) fClusterResNB[i] = recoParam->GetDefaultNonBendingReso(i);
 399     if (fClusterResB[i] < 0.) fClusterResB[i] = recoParam->GetDefaultBendingReso(i);
 400   }
 401
 402   // load geometry for track extrapolation to vertex and mono-cathod cluster finding
 403   if (!AliGeomManager::GetGeometry()) {
 404
 405     if (fReAlign) {
 406
 407       // get original geometry transformer
 408       AliGeomManager::LoadGeometry();
 409       if (!AliGeomManager::GetGeometry()) return;
 410       if (fOldAlignStorage != "none" && !AliGeomManager::ApplyAlignObjsFromCDB("MUON")) return;
 411       fOldGeoTransformer = new AliMUONGeometryTransformer();
 412       fOldGeoTransformer->LoadGeometryData();
 413
 414       // load the new geometry
 415       if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
 416       AliGeomManager::GetGeometry()->UnlockGeometry();
 417       AliGeomManager::LoadGeometry();
 418       if (!AliGeomManager::GetGeometry()) return;
 419       if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
 420       if (!fNewAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
 421       else {
 422         // recover default storage full name (raw:// cannot be used to set specific storage)
 423         TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
 424         if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
 425         else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
 426         cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
 427       }
 428       if (!AliGeomManager::ApplyAlignObjsFromCDB("MUON")) return;
 429
 430     } else {
 431
 432       AliGeomManager::LoadGeometry();
 433       if (!AliGeomManager::GetGeometry()) return;
 434       if (!AliGeomManager::ApplyAlignObjsFromCDB("MUON")) return;
 435
 436     }
 437
 438   } else fReAlign = kFALSE; // disable the realignment if the geometry was already loaded
 439
 440   fNewGeoTransformer = new AliMUONGeometryTransformer();
 441   fNewGeoTransformer->LoadGeometryData();
 442
 443   fTriggerCircuit = new AliMUONTriggerCircuit(fNewGeoTransformer);
 444
 445   fESDInterface = new AliMUONESDInterface();
 446
 447   fRefitter = new AliMUONRefitter(recoParam);
 448   fRefitter->Connect(fESDInterface);
 449
 450   AliMUONVDigitStore *digitStore = AliMUONESDInterface::NewDigitStore();
 451   fTrackHitPattern = new AliMUONTrackHitPattern(recoParam, *fNewGeoTransformer, *digitStore, 0x0);
 452   delete digitStore; // assume the digitStore is never used when we use AliMUONTrackHitPattern
 453 }
 454
 455 //________________________________________________________________________
 456 void AliAnalysisTaskMuonRefit::Terminate(Option_t *)
 457 {
 458 }
 459
 460 //________________________________________________________________________
 461 void AliAnalysisTaskMuonRefit::ModifyCluster(AliMUONVCluster& cl)
 462 {
 463   /// Reset the cluster resolution to the one given to the task and change
 464   /// the cluster position according to the new alignment parameters if required
 465
 466   Double_t gX,gY,gZ,lX,lY,lZ;
 467
 468   // change their resolution
 469   cl.SetErrXY(fClusterResNB[cl.GetChamberId()], fClusterResB[cl.GetChamberId()]);
 470
 471   // change their position
 472   if (fReAlign) {
 473     gX = cl.GetX();
 474     gY = cl.GetY();
 475     gZ = cl.GetZ();
 476     fOldGeoTransformer->Global2Local(cl.GetDetElemId(),gX,gY,gZ,lX,lY,lZ);
 477     fNewGeoTransformer->Local2Global(cl.GetDetElemId(),lX,lY,lZ,gX,gY,gZ);
 478     cl.SetXYZ(gX,gY,gZ);
 479   }
 480
 481   // "remove" mono-cathod clusters on stations 3-4-5 if required
 482   // (to be done after moving clusters to the new position)
 483   if (fRemoveMonoCathCl && cl.GetChamberId() > 3) {
 484     Bool_t hasBending, hasNonBending;
 485     if (fCheckAllPads) CheckPads(&cl, hasBending, hasNonBending);
 486     else CheckPadsBelow(&cl, hasBending, hasNonBending);
 487     if (!hasNonBending) cl.SetErrXY(10., cl.GetErrY());
 488     if (!hasBending) cl.SetErrXY(cl.GetErrX(), 10.);
 489   }
 490
 491 }
 492
 493 //________________________________________________________________________
 494 void AliAnalysisTaskMuonRefit::CheckPads(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
 495 {
 496   /// Check that this cluster contains pads on both cathods
 497
 498   // reset
 499   hasBending = kFALSE;
 500   hasNonBending = kFALSE;
 501
 502   // loop over digits contained in the cluster
 503   for (Int_t iDigit = 0; iDigit < cl->GetNDigits(); iDigit++) {
 504
 505     Int_t manuId = AliMUONVDigit::ManuId(cl->GetDigitId(iDigit));
 506
 507     // check the location of the manu the digit belongs to
 508     if (manuId > 0) {
 509       if (manuId & AliMpConstants::ManuMask(AliMp::kNonBendingPlane)) hasNonBending = kTRUE;
 510       else hasBending = kTRUE;
 511     }
 512
 513   }
 514
 515 }
 516
 517 //________________________________________________________________________
 518 void AliAnalysisTaskMuonRefit::CheckPadsBelow(AliMUONVCluster *cl, Bool_t &hasBending, Bool_t &hasNonBending) const
 519 {
 520   /// Check that this cluster contains pads on both cathods just under its position
 521
 522   // reset
 523   hasBending = kFALSE;
 524   hasNonBending = kFALSE;
 525
 526   // get the cathod corresponding to the bending/non-bending plane
 527   Int_t deId = cl->GetDetElemId();
 528   AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(deId, kFALSE);
 529   if (!de) return;
 530   AliMp::CathodType cath1 = de->GetCathodType(AliMp::kBendingPlane);
 531   AliMp::CathodType cath2 = de->GetCathodType(AliMp::kNonBendingPlane);
 532
 533   // get the corresponding segmentation
 534   const AliMpVSegmentation* seg1 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath1);
 535   const AliMpVSegmentation* seg2 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath2);
 536   if (!seg1 || !seg2) return;
 537
 538   // get local coordinate of the cluster
 539   Double_t lX,lY,lZ;
 540   Double_t gX = cl->GetX();
 541   Double_t gY = cl->GetY();
 542   Double_t gZ = cl->GetZ();
 543   fNewGeoTransformer->Global2Local(deId,gX,gY,gZ,lX,lY,lZ);
 544
 545   // find pads below the cluster
 546   AliMpPad pad1 = seg1->PadByPosition(lX, lY, kFALSE);
 547   AliMpPad pad2 = seg2->PadByPosition(lX, lY, kFALSE);
 548
 549   // build their ID if pads are valid
 550   UInt_t padId1 = (pad1.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad1.GetManuId(), pad1.GetManuChannel(), cath1) : 0;
 551   UInt_t padId2 = (pad2.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad2.GetManuId(), pad2.GetManuChannel(), cath2) : 0;
 552
 553   // check if the cluster contains these pads
 554   for (Int_t iDigit = 0; iDigit < cl->GetNDigits(); iDigit++) {
 555
 556     UInt_t digitId = cl->GetDigitId(iDigit);
 557
 558     if (digitId == padId1) {
 559
 560       hasBending = kTRUE;
 561       if (hasNonBending) break;
 562
 563     } else if (digitId == padId2) {
 564
 565       hasNonBending = kTRUE;
 566       if (hasBending) break;
 567
 568     }
 569
 570   }
 571
 572 }
 573
 574 //________________________________________________________________________
 575 Bool_t AliAnalysisTaskMuonRefit::SetMagField() const
 576 {
 577   // Dealing with the magnetic field map assuming nominal current
 578   // Construct the mag field map from the data in GRP and the given field map
 579   // Set the global mag field instance
 580
 581   if (TGeoGlobalMagField::Instance()->IsLocked()) delete TGeoGlobalMagField::Instance();
 582
 583   AliGRPManager grpMan;
 584   if (!grpMan.ReadGRPEntry()) {
 585     AliError("failed to load GRP Data from OCDB");
 586     return kFALSE;
 587   }
 588
 589   const AliGRPObject *grpData = grpMan.GetGRPData();
 590   if (!grpData) {
 591     AliError("GRP Data is not loaded");
 592     return kFALSE;
 593   }
 594
 595   Char_t l3Polarity = grpData->GetL3Polarity();
 596   if (l3Polarity == AliGRPObject::GetInvalidChar()) {
 597     AliError("GRP/GRP/Data entry:  missing value for the L3 polarity !");
 598     return kFALSE;
 599   }
 600
 601   Char_t diPolarity = grpData->GetDipolePolarity();
 602   if (diPolarity == AliGRPObject::GetInvalidChar()) {
 603     AliError("GRP/GRP/Data entry:  missing value for the dipole polarity !");
 604     return kFALSE;
 605   }
 606
 607   TString beamType = grpData->GetBeamType();
 608   if (beamType==AliGRPObject::GetInvalidString()) {
 609     AliError("GRP/GRP/Data entry:  missing value for the beam type !");
 610     return kFALSE;
 611   }
 612
 613   Float_t beamEnergy = grpData->GetBeamEnergy();
 614   if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
 615     AliError("GRP/GRP/Data entry:  missing value for the beam energy !");
 616     return kFALSE;
 617   }
 618
 619   AliMagF::BeamType_t btype = AliMagF::kNoBeamField;
 620   TString btypestr = beamType;
 621   btypestr.ToLower();
 622   TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1");
 623   TPRegexp ionBeam("(lead|pb|ion|a|A)\\s*-?\\s*\\1");
 624   TPRegexp protonionBeam("(proton|p)\\s*-?\\s*(lead|pb|ion|a|A)");
 625   TPRegexp ionprotonBeam("(lead|pb|ion|a|A)\\s*-?\\s*(proton|p)");
 626   if (btypestr.Contains(ionBeam)) btype = AliMagF::kBeamTypeAA;
 627   else if (btypestr.Contains(protonBeam)) btype = AliMagF::kBeamTypepp;
 628   else if (btypestr.Contains(protonionBeam)) btype = AliMagF::kBeamTypepA;
 629   else if (btypestr.Contains(ionprotonBeam)) btype = AliMagF::kBeamTypeAp;
 630   else AliInfoGeneral("AliMagF",Form("Assume no LHC magnet field for the beam type %s !",beamType.Data()));
 631
 632   TGeoGlobalMagField::Instance()->SetField(new AliMagF("fld", "fld", l3Polarity ? -1:1, diPolarity ? -1:1,
 633                                                        AliMagF::k5kG, btype, beamEnergy, 2, 15., fField.Data()));
 634   TGeoGlobalMagField::Instance()->Lock();
 635   AliInfo("Running with the B field constructed out of GRP and custom field map !");
 636
 637   return kTRUE;
 638
 639 }
 640