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