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 <TList.h>
  25 #include <TObjString.h>
  26 #include <TTree.h>
  27 #include <TString.h>
  28
  29 #include <AliCDBEntry.h>
  30 #include <AliCDBId.h>
  31 #include <AliCDBManager.h>
  32 #include <AliTRDCalDet.h>
  33
  34 #include <AliLog.h>
  35 #include <TTree.h>
  36 #include <TChain.h>
  37 #include <AliPID.h>
  38 #include <AliVEvent.h>
  39 #include <AliESDEvent.h>
  40 #include <AliESDpid.h>
  41 #include <AliESDtrack.h>
  42 #include <AliESDInputHandler.h>
  43 #include <AliAnalysisManager.h>
  44 #include <AliTRDPIDReference.h>
  45 #include <AliTRDPIDResponse.h>
  46 #include <AliTender.h>
  47
  48 #include "AliTRDTenderSupply.h"
  49
  50 ClassImp(AliTRDTenderSupply)
  51
  52 //_____________________________________________________
  53 AliTRDTenderSupply::AliTRDTenderSupply() :
  54   AliTenderSupply(),
  55   fESD(NULL),
  56   fESDpid(NULL),
  57   fChamberGainOld(NULL),
  58   fChamberGainNew(NULL),
  59   fChamberVdriftOld(NULL),
  60   fChamberVdriftNew(NULL),
  61   fRefFilename(""),
  62   fPIDmethod(kNNpid),
  63   fNormalizationFactor(1.),
  64   fPthreshold(0.8),
  65   fNBadChambers(0),
  66   fGainCorrection(kTRUE),
  67   fLoadReferencesFromCDB(kFALSE),
  68   fHasReferences(kFALSE)
  69 {
  70   //
  71   // default ctor
  72   //
  73 }
  74
  75 //_____________________________________________________
  76 AliTRDTenderSupply::AliTRDTenderSupply(const char *name, const AliTender *tender) :
  77   AliTenderSupply(name,tender),
  78   fESD(NULL),
  79   fESDpid(NULL),
  80   fChamberGainOld(NULL),
  81   fChamberGainNew(NULL),
  82   fChamberVdriftOld(NULL),
  83   fChamberVdriftNew(NULL),
  84   fRefFilename(""),
  85   fPIDmethod(kNNpid),
  86   fNormalizationFactor(1.),
  87   fPthreshold(0.8),
  88   fNBadChambers(0),
  89   fGainCorrection(kTRUE),
  90   fLoadReferencesFromCDB(kFALSE),
  91   fHasReferences(kFALSE)
  92 {
  93   //
  94   // named ctor
  95   //
  96   memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers);
  97 }
  98
  99 //_____________________________________________________
 100 AliTRDTenderSupply::~AliTRDTenderSupply()
 101 {
 102   //
 103   // dtor
 104   //
 105 }
 106
 107 //_____________________________________________________
 108 void AliTRDTenderSupply::Init()
 109 {
 110   //
 111   // Initialise TRD tender
 112   //
 113
 114   AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
 115
 116   // 1DLQ PID implemented in the AliESD object
 117   fESDpid=fTender->GetESDhandler()->GetESDpid();
 118   if (!fESDpid) {
 119     fESDpid=new AliESDpid;
 120     fTender->GetESDhandler()->SetESDpid(fESDpid);
 121   }
 122   // Load References
 123   if(!fLoadReferencesFromCDB) LoadReferences();
 124   fESDpid->GetTRDResponse().SetGainNormalisationFactor(fNormalizationFactor);
 125
 126   // Set Normalisation Factors
 127   if(mgr->GetMCtruthEventHandler()){
 128     // Assume MC
 129     //fESDpid->GetTRDResponse().SetGainNormalisationFactor(1.);
 130     SwitchOffGainCorrection();
 131   }
 132   else{
 133     // Assume Data
 134     //if(fPIDmethod == kNNpid) fPidRecal->SetGainScaleFactor(1.14);
 135     //fESDpid->GetTRDResponse().SetGainNormalisationFactor(1.14);
 136     SwitchOnGainCorrection();
 137   }
 138 }
 139
 140 //_____________________________________________________
 141 void AliTRDTenderSupply::ProcessEvent()
 142 {
 143   //
 144   // Reapply pid information
 145   //
 146   if (fTender->RunChanged()){
 147     AliDebug(0, Form("AliTPCTenderSupply::ProcessEvent - Run Changed (%d)\n",fTender->GetRun()));
 148     if (fGainCorrection) SetChamberGain();
 149     if(!fHasReferences) LoadReferences();
 150   }
 151
 152   fESD = fTender->GetEvent();
 153   if (!fESD) return;
 154   Int_t ntracks=fESD->GetNumberOfTracks();
 155
 156   //
 157   // recalculate PID probabilities
 158   //
 159   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 160     AliESDtrack *track=fESD->GetTrack(itrack);
 161     // Recalculate likelihoods
 162     if(!(track->GetStatus() & AliESDtrack::kTRDout)) continue;
 163     if(fGainCorrection) ApplyGainCorrection(track);
 164     switch(fPIDmethod){
 165       case kNNpid:
 166         break;
 167       case k1DLQpid:
 168         fESDpid->MakeTRDPID(fESD->GetTrack(itrack));
 169         break;
 170       default:
 171         AliError("PID Method not implemented (yet)");
 172     }
 173   }
 174 }
 175
 176 //_____________________________________________________
 177 void AliTRDTenderSupply::LoadReferences(){
 178   //
 179   // Load Reference from the OCDB/OADB into the PID Response
 180   //
 181   if(fLoadReferencesFromCDB){
 182     AliDebug(1, "Loading Reference Distributions from the OCDB");
 183     AliCDBEntry *en = fTender->GetCDBManager()->Get("TRD/Calib/PIDLQ1D");
 184     if(!en){
 185       AliError("References for 1D Likelihood Method not in OCDB");
 186       return;
 187     }
 188     en->GetId().Print();
 189     TObjArray *arr = dynamic_cast<TObjArray *>(en->GetObject());
 190     if(!arr) AliError("List with the references not found");
 191
 192     // Get new references
 193     TIter refs(arr);
 194     TObject *o = NULL;
 195     AliTRDPIDReference *ref = NULL;
 196     while((o = refs())){
 197       if(!TString(o->IsA()->GetName()).CompareTo("AliTRDPIDReference")){
 198         ref = dynamic_cast<AliTRDPIDReference *>(o);
 199         break;
 200       }
 201     }
 202     if(ref){
 203       fESDpid->GetTRDResponse().Load(ref);
 204       fHasReferences = kTRUE;
 205       AliInfo("Reference distributions loaded into the PID Response");
 206     } else {
 207       AliError("References not found");
 208     }
 209   } else {
 210     // Backward compatibility mode
 211     AliInfo("Loading Reference Distributions from ROOT file");
 212     if(fRefFilename.Length() != 0){
 213         fESDpid->GetTRDResponse().Load(fRefFilename.Data());
 214         fHasReferences = kTRUE;
 215     } else{
 216         AliError("No file defined");
 217     }
 218   }
 219 }
 220
 221 //_____________________________________________________
 222 void AliTRDTenderSupply::SetChamberGain(){
 223   //
 224   // Load Chamber Gain factors into the Tender supply
 225   //
 226
 227   //find previous entry from the UserInfo
 228   TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree();
 229   if (!tree) {
 230   AliError("Tree not found in ESDhandler");
 231     return;
 232   }
 233
 234   TList *userInfo=(TList*)tree->GetUserInfo();
 235   if (!userInfo) {
 236     AliError("No UserInfo found in tree");
 237     return;
 238   }
 239
 240   TList *cdbList=(TList*)userInfo->FindObject("cdbList");
 241   if (!cdbList) {
 242     AliError("No cdbList found in UserInfo");
 243     if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
 244     return;
 245   }
 246
 247   TIter nextCDB(cdbList);
 248   TObjString *os=0x0;
 249   while ( (os=(TObjString*)nextCDB()) ){
 250     if(os->GetString().Contains("TRD/Calib/ChamberGainFactor")){
 251       // Get Old gain calibration
 252       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
 253
 254       AliCDBEntry *entry=fTender->GetCDBManager()->Get(id->GetPath(), id->GetFirstRun(), id->GetVersion());
 255       if (!entry) {
 256         AliError("No previous gain calibration entry found");
 257         return;
 258       }
 259
 260       fChamberGainOld = dynamic_cast<AliTRDCalDet *>(entry->GetObject());
 261
 262       AliDebug(1, Form("Used old Gain entry: %s\n",entry->GetId().ToString().Data()));
 263     } else if(os->GetString().Contains("TRD/Calib/ChamberVdrift")){
 264       // Get Old drift velocity calibration
 265       AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
 266
 267       AliCDBEntry *entry=fTender->GetCDBManager()->Get(id->GetPath(), id->GetFirstRun(), id->GetVersion());
 268       if (!entry) {
 269         AliError("No previous drift velocity calibration entry found");
 270         return;
 271       }
 272
 273       fChamberVdriftOld = dynamic_cast<AliTRDCalDet *>(entry->GetObject());
 274
 275       AliDebug(1, Form("Used old Drift velocity entry: %s\n",entry->GetId().ToString().Data()));
 276
 277     }
 278   }
 279
 280   // Get Latest Gain Calib Object
 281   AliCDBEntry *entryNew=fTender->GetCDBManager()->Get("TRD/Calib/ChamberGainFactor",fTender->GetRun());
 282   if (entryNew) {
 283     AliDebug(1, Form("Used new Gain entry: %s\n",entryNew->GetId().ToString().Data()));
 284     fChamberGainNew = dynamic_cast<AliTRDCalDet *>(entryNew->GetObject());
 285   } else
 286     AliError("No new gain calibration entry found");
 287
 288   // Also get the latest Drift Velocity calibration object
 289   entryNew=fTender->GetCDBManager()->Get("TRD/Calib/ChamberVdrift",fTender->GetRun());
 290   if (entryNew) {
 291     AliDebug(1, Form("Used new Drift velocity entry: %s\n",entryNew->GetId().ToString().Data()));
 292     fChamberVdriftNew = dynamic_cast<AliTRDCalDet *>(entryNew->GetObject());
 293   } else
 294     AliError("No new drift velocity calibration entry found");
 295
 296   if(!fChamberGainNew || !fChamberVdriftNew) AliError("No recent calibration found");
 297 }
 298
 299 //_____________________________________________________
 300 void AliTRDTenderSupply::ApplyGainCorrection(AliESDtrack * track){
 301   //
 302   // Apply new gain factors to the track
 303   //
 304   if(!fChamberGainNew || !fChamberGainOld){
 305     AliError("Cannot apply gain correction.");
 306     return;
 307   }
 308
 309   if(!(track->GetStatus() & AliESDtrack::kTRDout)) return;
 310   Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
 311   if(p < fPthreshold) return; // Apply low momentum cutoff
 312
 313   Bool_t applyCorrectionVdrift = kFALSE;
 314   if(fChamberVdriftOld && fChamberVdriftNew) applyCorrectionVdrift = kTRUE;
 315
 316   Int_t chamberID[kNPlanes];
 317   for(Int_t il = 0; il < kNPlanes; il++) chamberID[il] = -1;
 318   if(!GetTRDchamberID(track, chamberID)) return;
 319   Int_t nTrackletsPID = 0, nTracklets = track->GetTRDntracklets();
 320   for(Int_t iplane = 0; iplane < kNPlanes; iplane++){
 321     if(chamberID[iplane] < 0) continue;
 322
 323     // Mask out bad chambers
 324     Bool_t isMasked = kFALSE;
 325     for(UInt_t icam = 0; icam < fNBadChambers; icam++)
 326       if(fBadChamberID[icam] == chamberID[iplane]){
 327         isMasked = kTRUE;
 328         break;
 329       }
 330     if(isMasked){
 331       for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
 332         track->SetTRDslice(0, iplane, islice);
 333       }
 334       continue;
 335     }
 336
 337     // Take old and new gain factor and make ratio
 338     Double_t facOld = fChamberGainOld->GetValue(chamberID[iplane]);
 339     Double_t facNew = fChamberGainNew->GetValue(chamberID[iplane]);
 340     Double_t correction = facNew/facOld;
 341     if(applyCorrectionVdrift){
 342       // apply also correction for drift velocity calibration
 343       Double_t vDriftOld = fChamberVdriftOld->GetValue(chamberID[iplane]);
 344       Double_t vDriftNew = fChamberVdriftNew->GetValue(chamberID[iplane]);
 345       correction *= vDriftNew/vDriftOld;
 346     }
 347     AliDebug(2, Form("Applying correction factor %f\n", correction));
 348     Int_t nSlices = 0;
 349     for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
 350       Double_t qslice = track->GetTRDslice(iplane, islice);
 351       if(qslice <= 0.) continue;
 352       track->SetTRDslice(qslice / correction, iplane, islice);
 353       nSlices++;
 354     }
 355     if(nSlices) nTrackletsPID++;
 356   }
 357   // Use nTrackletsPID to indicate the number of tracklets from good
 358   // chambers so they are used for the PID
 359   track->SetTRDntracklets(nTrackletsPID | (nTracklets << 3));
 360 }
 361
 362 //_____________________________________________________
 363 Bool_t AliTRDTenderSupply::GetTRDchamberID(AliESDtrack * const track, Int_t *detectors) {
 364   //
 365   // Calculate TRD chamber ID
 366   //
 367   Double_t xLayer[kNPlanes] = {300.2, 312.8, 325.4, 338., 350.6, 363.2};
 368   Double_t etamin[kNStacks] = {0.536, 0.157, -0.145, -0.527,-0.851};
 369   Double_t etamax[kNStacks] = {0.851, 0.527, 0.145, -0.157,-0.536};
 370   for(Int_t ily = 0; ily < kNPlanes; ily++) detectors[ily] = -1;
 371
 372   const AliExternalTrackParam *trueparam = NULL;
 373   if(track->GetTPCInnerParam()) trueparam = track->GetTPCInnerParam();
 374   else if(track->GetOuterParam()) trueparam = track->GetOuterParam();
 375   else if(track->GetInnerParam()) trueparam = track->GetInnerParam();
 376   if(!trueparam){
 377     AliDebug(2, "No Track Params");
 378     return kFALSE;
 379   }
 380
 381   AliExternalTrackParam workparam(*trueparam); // Do calculation on working Copy
 382   Double_t pos[3];
 383   Int_t nDet = 0;
 384   for(Int_t ily = 0; ily < kNPlanes; ily++){
 385     if(!workparam.PropagateTo(xLayer[ily], fESD->GetMagneticField())) {
 386       AliDebug(2, "Propagation failed");
 387       break;
 388     }
 389     workparam.GetXYZ(pos);
 390     Double_t trackAlpha = TMath::ATan2(pos[1], pos[0]);
 391     if(trackAlpha < 0) trackAlpha = 2 * TMath::Pi() + trackAlpha;
 392     Double_t secAlpha = 2 * TMath::Pi() / 18.;
 393
 394     Int_t sector = static_cast<Int_t>(trackAlpha/secAlpha);
 395     Double_t etaTrack = track->Eta();
 396     Int_t stack = -1;
 397     for(Int_t istack = 0; istack < 5; istack++){
 398       if(etaTrack >= etamin[istack] && etaTrack <= etamax[istack]){
 399         stack = istack;
 400         break;
 401       }
 402     }
 403     if(stack < 0) {
 404       AliDebug(2, "Dead Area");
 405       continue;
 406     }
 407
 408     detectors[ily] = sector * kNStacks * kNPlanes + stack * kNPlanes + ily;
 409     nDet++;
 410   }
 411   return nDet ? kTRUE : kFALSE;
 412 }
 413