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