ANALYSIS/TenderSupplies/AliTRDTenderSupply.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
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 // TRD tender: reapply pid on the fly                                        //
  20 //                                                                           //
  21 ///////////////////////////////////////////////////////////////////////////////
  22
  23 #include <TChain.h>
  24 #include <TDirectory.h>
  25 #include <TFile.h>
  26 #include <TList.h>
  27 #include <TObjString.h>
  28 #include <TTree.h>
  29 #include <TString.h>
  30 #include <TVectorD.h>
  31
  32 #include <AliCDBEntry.h>
  33 #include <AliCDBId.h>
  34 #include <AliCDBManager.h>
  35 #include <AliOADBContainer.h>
  36 #include <AliTRDCalDet.h>
  37
  38 #include <AliLog.h>
  39 #include <TTree.h>
  40 #include <TChain.h>
  41 #include <AliGeomManager.h>
  42 #include <AliPID.h>
  43 #include <AliVEvent.h>
  44 #include <AliESDEvent.h>
  45 #include <AliESDpid.h>
  46 #include <AliESDtrack.h>
  47 #include <AliESDInputHandler.h>
  48 #include <AliAnalysisManager.h>
  49 #include <AliTrackerBase.h>
  50 #include <AliTRDPIDResponseObject.h>
  51 #include <AliTRDPIDResponse.h>
  52 #include "AliTRDCalChamberStatus.h"
  53 #include <AliTender.h>
  54
  55 #include "AliTRDTenderSupply.h"
  56
  57 ClassImp(AliTRDTenderSupply)
  58
  59 //_____________________________________________________
  60 AliTRDTenderSupply::AliTRDTenderSupply() :
  61   AliTenderSupply(),
  62   fESD(NULL),
  63   fESDpid(NULL),
  64   fChamberGainOld(NULL),
  65   fChamberGainNew(NULL),
  66   fChamberVdriftOld(NULL),
  67   fChamberVdriftNew(NULL),
  68   fRunByRunCorrection(NULL),
  69   fPIDmethod(kNNpid),
  70   fNormalizationFactor(1.),
  71   fPthreshold(0.8),
  72   fNBadChambers(0),
  73   fGeoFile(NULL),
  74   fGainCorrection(kTRUE),
  75   fLoadReferences(kFALSE),
  76   fLoadReferencesFromCDB(kFALSE),
  77   fLoadDeadChambers(kFALSE),
  78   fHasReferences(kFALSE),
  79   fHasNewCalibration(kTRUE),
  80   fDebugMode(kFALSE),
  81   fNameRunByRunCorrection()
  82 {
  83   //
  84   // default ctor
  85   //
  86   memset(fSlicesForPID, 0, sizeof(UInt_t) * 2);
  87   memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers);
  88 }
  89
  90 //_____________________________________________________
  91 AliTRDTenderSupply::AliTRDTenderSupply(const char *name, const AliTender *tender) :
  92   AliTenderSupply(name,tender),
  93   fESD(NULL),
  94   fESDpid(NULL),
  95   fChamberGainOld(NULL),
  96   fChamberGainNew(NULL),
  97   fChamberVdriftOld(NULL),
  98   fChamberVdriftNew(NULL),
  99   fRunByRunCorrection(NULL),
 100   fPIDmethod(kNNpid),
 101   fNormalizationFactor(1.),
 102   fPthreshold(0.8),
 103   fNBadChambers(0),
 104   fGeoFile(NULL),
 105   fGainCorrection(kTRUE),
 106   fLoadReferences(kFALSE),
 107   fLoadReferencesFromCDB(kFALSE),
 108   fLoadDeadChambers(kFALSE),
 109   fHasReferences(kFALSE),
 110   fHasNewCalibration(kTRUE),
 111   fDebugMode(kFALSE),
 112   fNameRunByRunCorrection()
 113 {
 114   //
 115   // named ctor
 116   //
 117   memset(fSlicesForPID, 0, sizeof(UInt_t) * 2);
 118   memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers);
 119 }
 120
 121 //_____________________________________________________
 122 AliTRDTenderSupply::~AliTRDTenderSupply()
 123 {
 124   //
 125   // dtor
 126   //
 127 }
 128
 129 //_____________________________________________________
 130 void AliTRDTenderSupply::Init()
 131 {
 132   //
 133   // Initialise TRD tender
 134   //
 135
 136   AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
 137
 138   // 1DLQ PID implemented in the AliESD object
 139   fESDpid=fTender->GetESDhandler()->GetESDpid();
 140   if (!fESDpid) {
 141     fESDpid=new AliESDpid;
 142     fTender->GetESDhandler()->SetESDpid(fESDpid);
 143   }
 144   // Load References
 145   if(fLoadReferences && !fLoadReferencesFromCDB) LoadReferences();
 146   //fESDpid->GetTRDResponse().SetGainNormalisationFactor(fNormalizationFactor);
 147   fESDpid->SetTRDslicesForPID(fSlicesForPID[0], fSlicesForPID[1]);
 148
 149   if(fNameRunByRunCorrection.Length()) LoadRunByRunCorrection(fNameRunByRunCorrection.Data());
 150
 151   // Set Normalisation Factors
 152   if(mgr->GetMCtruthEventHandler()){
 153     // Assume MC
 154     //fESDpid->GetTRDResponse().SetGainNormalisationFactor(1.);
 155     SwitchOffGainCorrection();
 156   }
 157   else{
 158     // Assume Data
 159     //if(fPIDmethod == kNNpid) fPidRecal->SetGainScaleFactor(1.14);
 160     //fESDpid->GetTRDResponse().SetGainNormalisationFactor(1.14);
 161     SwitchOnGainCorrection();
 162   }
 163 }
 164
 165 //_____________________________________________________
 166 void AliTRDTenderSupply::ProcessEvent()
 167 {
 168   //
 169   // Reapply pid information
 170   //
 171   if (fTender->RunChanged()){
 172     AliDebug(0, Form("AliTPCTenderSupply::ProcessEvent - Run Changed (%d)\n",fTender->GetRun()));
 173     if (fGainCorrection) SetChamberGain();
 174     if(fLoadReferences && !fHasReferences) LoadReferences();
 175     if(fLoadDeadChambers) LoadDeadChambersFromCDB();
 176     // Load Geometry
 177     if(AliGeomManager::GetGeometry()){
 178       AliInfo("Geometry already loaded by other tenders");
 179     } else {
 180       if(fGeoFile) AliInfo(Form("Load geometry from file %s\n", fGeoFile));
 181       else AliInfo("Load Geometry from OCDB\n");
 182       AliGeomManager::LoadGeometry(fGeoFile);
 183     }
 184   }
 185
 186
 187   fESD = fTender->GetEvent();
 188   if (!fESD) return;
 189   Int_t ntracks=fESD->GetNumberOfTracks();
 190
 191   //
 192   // recalculate PID probabilities
 193   //
 194   Int_t detectors[kNPlanes];
 195   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 196     for(Int_t idet = 0; idet < 5; idet++) detectors[idet] = -1;
 197     AliESDtrack *track=fESD->GetTrack(itrack);
 198     // Recalculate likelihoods
 199     if(!(track->GetStatus() & AliESDtrack::kTRDout)) continue;
 200     AliDebug(2, Form("TRD track found, gain correction: %s, Number of bad chambers: %d\n", fGainCorrection ? "Yes" : "No", fNBadChambers));
 201     if(GetTRDchamberID(track, detectors)){
 202       if(fGainCorrection && fHasNewCalibration) ApplyGainCorrection(track, detectors);
 203       if(fNBadChambers) MaskChambers(track, detectors);
 204     }
 205     if(fRunByRunCorrection) ApplyRunByRunCorrection(track);
 206     if(fNormalizationFactor != 1.){
 207       //printf("Gain Factor: %f\n", fNormalizationFactor);
 208       // Renormalize charge
 209       Double_t qslice = -1;
 210       for(Int_t ily = 0; ily < 6; ily++){
 211         for(Int_t is = 0; is < track->GetNumberOfTRDslices(); is++){
 212           qslice = track->GetTRDslice(ily, is);
 213           //printf("Doing layer %d slice %d, value %f\n", ily, is, qslice);
 214           if(qslice >0){
 215             qslice *= fNormalizationFactor;
 216             //printf("qslice new: %f\n", qslice);
 217             track->SetTRDslice(qslice, ily, is);
 218           }
 219         }
 220       }
 221     }
 222     switch(fPIDmethod){
 223       case kNNpid:
 224         break;
 225       case k1DLQpid:
 226         fESDpid->MakeTRDPID(fESD->GetTrack(itrack));
 227         break;
 228       default:
 229         AliError("PID Method not implemented (yet)");
 230     }
 231   }
 232 }
 233
 234 //_____________________________________________________
 235 void AliTRDTenderSupply::LoadDeadChambersFromCDB(){
 236   //
 237   // Load Dead Chambers from the OCDB
 238   //
 239   AliDebug(1, "Loading Dead Chambers from the OCDB");
 240   AliCDBEntry *en = fTender->GetCDBManager()->Get("TRD/Calib/ChamberStatus",fTender->GetRun());
 241   if(!en){
 242    AliError("Dead Chambers not in OCDB");
 243    return;
 244   }
 245   en->GetId().Print();
 246
 247   AliTRDCalChamberStatus* chamberStatus = 0;
 248   if(en){
 249     chamberStatus = (AliTRDCalChamberStatus*)en->GetObject();
 250     if(!chamberStatus) AliError("List with the dead chambers not found");
 251     for(Int_t ichamber = 0; ichamber < 540; ichamber++) {
 252       if(!chamberStatus->IsGood(ichamber)){
 253         //printf("Chamber not installed %d\n",ichamber);
 254         AddBadChamber(ichamber);
 255       }
 256     }
 257   }
 258 }
 259
 260 //_____________________________________________________
 261 void AliTRDTenderSupply::LoadReferences(){
 262   //
 263   // Load Reference from the OCDB/OADB into the PID Response
 264   //
 265   if(fLoadReferencesFromCDB){
 266     AliDebug(1, "Loading Reference Distributions from the OCDB");
 267     AliCDBEntry *en = fTender->GetCDBManager()->Get("TRD/Calib/PIDLQ1D");
 268     if(!en){
 269       AliError("References for 1D Likelihood Method not in OCDB");
 270       return;
 271     }
 272     en->GetId().Print();
 273     TObjArray *arr = dynamic_cast<TObjArray *>(en->GetObject());
 274     if(!arr) AliError("List with the references not found");
 275
 276     // Get new references
 277     TIter refs(arr);
 278     TObject *o = NULL;
 279     AliTRDPIDResponseObject *ref = NULL;
 280     while((o = refs())){
 281       if(!TString(o->IsA()->GetName()).CompareTo("AliTRDPIDResponseObject")){
 282         ref = dynamic_cast<AliTRDPIDResponseObject *>(o);
 283         break;
 284       }
 285     }
 286     if(ref){
 287       fESDpid->GetTRDResponse().SetPIDResponseObject(ref);
 288       fHasReferences = kTRUE;
 289       AliDebug(1, "Reference distributions loaded into the PID Response");
 290     } else {
 291       AliError("References not found");
 292     }
 293   } else {
 294     // Backward compatibility mode
 295     AliInfo("Loading Reference Distributions from ROOT file");
 296     fESDpid->GetTRDResponse().Load("$TRAIN_ROOT/util/tender/LQ1dRef_v3.root");
 297     fHasReferences = kTRUE;
 298   }
 299 }
 300
 301 //_____________________________________________________
 302 void AliTRDTenderSupply::SetChamberGain(){
 303   //
 304   // Load Chamber Gain factors into the Tender supply
 305   //
 306
 307   //find previous entry from the UserInfo
 308   TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree();
 309   if (!tree) {
 310     fHasNewCalibration = kFALSE;
 311     AliError("Tree not found in ESDhandler");
 312     return;
 313   }
 314
 315   TList *userInfo=(TList*)tree->GetUserInfo();
 316   if (!userInfo) {
 317     fHasNewCalibration = kFALSE;
 318     AliError("No UserInfo found in tree");
 319     return;
 320   }
 321
 322   TList *cdbList=(TList*)userInfo->FindObject("cdbList");
 323   if (!cdbList) {
 324     fHasNewCalibration = kFALSE;
 325     AliError("No cdbList found in UserInfo");
 326     if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
 327     return;
 328   }
 329   fHasNewCalibration = kTRUE;
 330
 331   TIter nextCDB(cdbList);
 332   TObjString *os=0x0;
 333   while ( (os=(TObjString*)nextCDB()) ){
 334     if(os->GetString().Contains("TRD/Calib/ChamberGainFactor")){
 335       // Get Old gain calibration
 336       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
 337
 338       AliCDBEntry *entry=fTender->GetCDBManager()->Get(id->GetPath(), id->GetFirstRun(), id->GetVersion());
 339       if (!entry) {
 340         AliError("No previous gain calibration entry found");
 341         return;
 342       }
 343
 344       fChamberGainOld = dynamic_cast<AliTRDCalDet *>(entry->GetObject());
 345
 346       AliDebug(1, Form("Used old Gain entry: %s\n",entry->GetId().ToString().Data()));
 347     } else if(os->GetString().Contains("TRD/Calib/ChamberVdrift")){
 348       // Get Old drift velocity calibration
 349       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
 350
 351       AliCDBEntry *entry=fTender->GetCDBManager()->Get(id->GetPath(), id->GetFirstRun(), id->GetVersion());
 352       if (!entry) {
 353         AliError("No previous drift velocity calibration entry found");
 354         return;
 355       }
 356
 357       fChamberVdriftOld = dynamic_cast<AliTRDCalDet *>(entry->GetObject());
 358
 359       AliDebug(1, Form("Used old Drift velocity entry: %s\n",entry->GetId().ToString().Data()));
 360
 361     }
 362   }
 363
 364   // Get Latest Gain Calib Object
 365   AliCDBEntry *entryNew=fTender->GetCDBManager()->Get("TRD/Calib/ChamberGainFactor",fTender->GetRun());
 366   if (entryNew) {
 367     AliDebug(1, Form("Used new Gain entry: %s\n",entryNew->GetId().ToString().Data()));
 368     fChamberGainNew = dynamic_cast<AliTRDCalDet *>(entryNew->GetObject());
 369   } else
 370     AliError("No new gain calibration entry found");
 371
 372   // Also get the latest Drift Velocity calibration object
 373   entryNew=fTender->GetCDBManager()->Get("TRD/Calib/ChamberVdrift",fTender->GetRun());
 374   if (entryNew) {
 375     AliDebug(1, Form("Used new Drift velocity entry: %s\n",entryNew->GetId().ToString().Data()));
 376     fChamberVdriftNew = dynamic_cast<AliTRDCalDet *>(entryNew->GetObject());
 377   } else
 378     AliError("No new drift velocity calibration entry found");
 379
 380   if(!fChamberGainNew || !fChamberVdriftNew){
 381     AliError("No recent calibration found");
 382     fHasNewCalibration = kFALSE;
 383   }
 384 }
 385
 386 //_____________________________________________________
 387 void AliTRDTenderSupply::LoadRunByRunCorrection(const char *filename){
 388   //
 389   // Define run by run gain correction for the charge
 390   //
 391
 392   TDirectory *bkp = gDirectory;
 393   TFile *in = TFile::Open(filename);
 394   bkp->cd();
 395   fRunByRunCorrection = dynamic_cast<AliOADBContainer *>(in->Get("TRDchargeCorrection"));
 396   delete in;
 397   if(fRunByRunCorrection )
 398     AliDebug(2, Form("OADB Container has %d runs\n", fRunByRunCorrection->GetNumberOfEntries()));
 399   /* Temporarily out due to a bug in AliOADBContainer
 400   fRunByRunCorrection = new AliOADBContainer("TRDchargeCorrection");
 401   Int_t status = fRunByRunCorrection->InitFromFile(filename, "TRDchargeCorrection");
 402   if(!status) AliDebug(1, Form("Run-dependend gain correction factors loaded from OADB file %s", filename));
 403   else{
 404     AliDebug(1, "Failed Loading Run-dependend gain correction factors");
 405     delete fRunByRunCorrection;
 406     fRunByRunCorrection = NULL;
 407   }
 408   */
 409 }
 410
 411 //_____________________________________________________
 412 Bool_t AliTRDTenderSupply::IsBadChamber(Int_t chamberID){
 413   //
 414   // Check if the chamber id is in the list of bad chambers
 415   //
 416   Bool_t isBad = kFALSE;
 417   for(UInt_t icam = 0; icam < fNBadChambers; icam++)
 418     if(fBadChamberID[icam] == chamberID){
 419             isBad = kTRUE;
 420             //printf("cross checking: %i \n",chamberID);
 421       break;
 422     }
 423   return isBad;
 424 }
 425
 426 //_____________________________________________________
 427 void AliTRDTenderSupply::ApplyGainCorrection(AliESDtrack * track, const Int_t * const chamberID){
 428   //
 429   // Apply new gain factors to the track
 430   //
 431   if(!fChamberGainNew || !fChamberGainOld){
 432     AliError("Cannot apply gain correction.");
 433     return;
 434   }
 435
 436   if(!(track->GetStatus() & AliESDtrack::kTRDout)) return;
 437   Bool_t applyCorrectionVdrift = kFALSE;
 438   if(fChamberVdriftOld && fChamberVdriftNew) applyCorrectionVdrift = kTRUE;
 439
 440   for(Int_t iplane = 0; iplane < kNPlanes; iplane++){
 441     if(chamberID[iplane] < 0) continue;
 442     if(IsBadChamber(chamberID[iplane])) continue; // Don't apply gain correction for chambers which are in the list of bad chambers
 443
 444     // Take old and new gain factor and make ratio
 445     Double_t facOld = fChamberGainOld->GetValue(chamberID[iplane]);
 446     Double_t facNew = fChamberGainNew->GetValue(chamberID[iplane]);
 447     Double_t correction = facNew/facOld;
 448     if(applyCorrectionVdrift){
 449       // apply also correction for drift velocity calibration
 450       Double_t vDriftOld = fChamberVdriftOld->GetValue(chamberID[iplane]);
 451       Double_t vDriftNew = fChamberVdriftNew->GetValue(chamberID[iplane]);
 452       correction *= vDriftNew/vDriftOld;
 453     }
 454     AliDebug(2, Form("Applying correction factor %f\n", correction));
 455     for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
 456       Double_t qslice = track->GetTRDslice(iplane, islice);
 457       if(qslice <= 0.) continue;
 458       track->SetTRDslice(qslice / correction, iplane, islice);
 459     }
 460   }
 461 }
 462
 463 //_____________________________________________________
 464 void AliTRDTenderSupply::ApplyRunByRunCorrection(AliESDtrack *const track) {
 465   //
 466   // Equalize charge distribution by applying run-by-run correction (multiplicative)
 467   //
 468
 469   TVectorD *corrfactor = dynamic_cast<TVectorD *>(fRunByRunCorrection->GetObject(fTender->GetRun()));
 470   if(!corrfactor) {
 471     AliDebug(2, "Couldn't derive gain correction factor from OADB");
 472     return;
 473   }
 474   else AliDebug(2, Form("Gain factor from OADB %f", (*corrfactor)[0]));
 475   Double_t slice = 0;
 476   for(Int_t ily = 0; ily < kNPlanes; ily++){
 477     for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
 478       slice = track->GetTRDslice(ily, islice);
 479       if(slice < 0.001) continue;   // do not modify slices which are 0 or negative
 480       slice *= (*corrfactor)[0];
 481       track->SetTRDslice(slice, ily, islice);
 482     }
 483   }
 484 }
 485
 486 //_____________________________________________________
 487 void AliTRDTenderSupply::MaskChambers(AliESDtrack *const track, const Int_t * const chamberID){
 488   //
 489   // Mask out chambers which are in the list of bad chambers
 490   // Set chamber signal to 0 and reduce the number of tracklets used for PID
 491   //
 492   AliDebug(2, "Masking bad chambers for TRD track");
 493   Int_t nTrackletsPID = 0, nslice = 0, nTracklets = track->GetTRDntracklets();
 494   Bool_t badChamber = kFALSE;
 495   //Int_t nbad = 0 ;      // Number of bad chambers which contain also a signal
 496   //Int_t nsliceBad = 0;  // Number of slices in tracklet in a bad chamber
 497   for(Int_t iplane = 0; iplane < kNPlanes; iplane++){
 498     badChamber = kFALSE;
 499     nslice = 0; //nsliceBad = 0;
 500     if(IsBadChamber(chamberID[iplane])) badChamber = kTRUE;
 501     for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
 502       if(badChamber){
 503         //if(track->GetTRDslice(iplane, islice)) nsliceBad++;
 504         track->SetTRDslice(-1, iplane, islice);
 505       } else if(track->GetTRDslice(iplane, islice) > 0.001) nslice++;
 506     }
 507     //if(nsliceBad) nbad++;
 508     if(nslice > 0) nTrackletsPID++;
 509   }
 510   //if(nbad) track->SetTRDncls(track->GetTRDncls() - 20 * nbad);      // subtract mean number of clusters per tracklet for bad tracklets
 511   // Use nTrackletsPID to indicate the number of tracklets from good
 512   // chambers so they are used for the PID
 513   track->SetTRDntracklets(nTrackletsPID | (nTracklets << 3));
 514 }
 515
 516 //_____________________________________________________
 517 Bool_t AliTRDTenderSupply::GetTRDchamberID(AliESDtrack * const track, Int_t *detectors) {
 518   //
 519   // Calculate TRD chamber ID
 520   //
 521   Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
 522   if(p < fPthreshold) return kFALSE; // Apply low momentum cutoff
 523
 524   Double_t xLayer[kNPlanes] = {300.2, 312.8, 325.4, 338., 350.6, 363.2};
 525   Double_t etamin[kNStacks] = {0.536, 0.157, -0.145, -0.527,-0.851};
 526   Double_t etamax[kNStacks] = {0.851, 0.527, 0.145, -0.157,-0.536};
 527   for(Int_t ily = 0; ily < kNPlanes; ily++) detectors[ily] = -1;
 528
 529   const AliExternalTrackParam *trueparam = NULL;
 530   if(track->GetOuterParam()) trueparam = track->GetOuterParam();
 531   else if(track->GetTPCInnerParam()) trueparam = track->GetTPCInnerParam();
 532   else if(track->GetInnerParam()) trueparam = track->GetInnerParam();
 533   if(!trueparam){
 534     AliDebug(2, "No Track Params");
 535     return kFALSE;
 536   }
 537
 538   AliExternalTrackParam workparam(*trueparam); // Do calculation on working Copy
 539   Double_t pos[3];
 540   Int_t nDet = 0;
 541   for(Int_t ily = 0; ily < kNPlanes; ily++){
 542     if(!AliTrackerBase::PropagateTrackToBxByBz(&workparam, xLayer[ily], 0.139, 100)){   // Assuming the pion mass
 543       AliDebug(2, "Propagation failed");
 544       break;
 545     }
 546     workparam.GetXYZ(pos);
 547     Double_t trackAlpha = TMath::ATan2(pos[1], pos[0]);
 548     if(fDebugMode){
 549       // Compare to simple propagation without magnetic field
 550       AliExternalTrackParam workparam1(*trueparam); // Do calculation on working Copy
 551       Double_t pos1[3];
 552       if(!workparam1.PropagateTo(xLayer[ily], fESD->GetMagneticField())) {
 553         AliDebug(2, "Propagation failed");
 554         break;
 555       }
 556       workparam1.GetXYZ(pos1);
 557       Double_t trackAlpha1 = TMath::ATan2(pos1[1], pos1[0]);
 558       AliDebug(2, Form("Alpha: Old %f, New %f, diff %f", trackAlpha1, trackAlpha, trackAlpha-trackAlpha1));
 559     }
 560     if(trackAlpha < 0) trackAlpha = 2 * TMath::Pi() + trackAlpha;
 561     Double_t secAlpha = 2 * TMath::Pi() / 18.;
 562
 563     Int_t sector = static_cast<Int_t>(trackAlpha/secAlpha);
 564     Double_t etaTrack = track->Eta();
 565     Int_t stack = -1;
 566     for(Int_t istack = 0; istack < 5; istack++){
 567       if(etaTrack >= etamin[istack] && etaTrack <= etamax[istack]){
 568         stack = istack;
 569         break;
 570       }
 571     }
 572     if(stack < 0) {
 573       AliDebug(2, "Dead Area");
 574       continue;
 575     }
 576
 577     detectors[ily] = sector * kNStacks * kNPlanes + stack * kNPlanes + ily;
 578     nDet++;
 579   }
 580   return nDet ? kTRUE : kFALSE;
 581 }
 582