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