Update of the TRD PID Response:
[u/mrichter/AliRoot.git] / ANALYSIS / TenderSupplies / AliTRDTenderSupply.cxx
CommitLineData
e75408ba 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
00685073 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>
e75408ba 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>
51a0ce25 44#include <AliTRDPIDReference.h>
e75408ba 45#include <AliTRDPIDResponse.h>
46#include <AliTender.h>
47
48#include "AliTRDTenderSupply.h"
49
00685073 50ClassImp(AliTRDTenderSupply)
e75408ba 51
00685073 52//_____________________________________________________
e75408ba 53AliTRDTenderSupply::AliTRDTenderSupply() :
54 AliTenderSupply(),
00685073 55 fESD(NULL),
56 fESDpid(NULL),
57 fChamberGainOld(NULL),
58 fChamberGainNew(NULL),
59 fChamberVdriftOld(NULL),
60 fChamberVdriftNew(NULL),
51a0ce25 61 fRefFilename(""),
9daf22f8 62 fPIDmethod(kNNpid),
51a0ce25 63 fNormalizationFactor(1.),
00685073 64 fPthreshold(0.8),
65 fNBadChambers(0),
51a0ce25 66 fGainCorrection(kTRUE),
67 fLoadReferencesFromCDB(kFALSE),
68 fHasReferences(kFALSE)
e75408ba 69{
70 //
71 // default ctor
72 //
73}
74
75//_____________________________________________________
76AliTRDTenderSupply::AliTRDTenderSupply(const char *name, const AliTender *tender) :
77 AliTenderSupply(name,tender),
00685073 78 fESD(NULL),
79 fESDpid(NULL),
80 fChamberGainOld(NULL),
81 fChamberGainNew(NULL),
82 fChamberVdriftOld(NULL),
83 fChamberVdriftNew(NULL),
51a0ce25 84 fRefFilename(""),
9daf22f8 85 fPIDmethod(kNNpid),
51a0ce25 86 fNormalizationFactor(1.),
00685073 87 fPthreshold(0.8),
88 fNBadChambers(0),
51a0ce25 89 fGainCorrection(kTRUE),
90 fLoadReferencesFromCDB(kFALSE),
91 fHasReferences(kFALSE)
e75408ba 92{
93 //
94 // named ctor
95 //
00685073 96 memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers);
e75408ba 97}
98
00685073 99//_____________________________________________________
e75408ba 100AliTRDTenderSupply::~AliTRDTenderSupply()
101{
102 //
103 // dtor
104 //
105}
106
107//_____________________________________________________
108void AliTRDTenderSupply::Init()
109{
110 //
111 // Initialise TRD tender
112 //
113
e75408ba 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 }
51a0ce25 122 // Load References
123 if(!fLoadReferencesFromCDB) LoadReferences();
124 fESDpid->GetTRDResponse().SetGainNormalisationFactor(fNormalizationFactor);
125
e75408ba 126 // Set Normalisation Factors
127 if(mgr->GetMCtruthEventHandler()){
128 // Assume MC
9daf22f8 129 //fESDpid->GetTRDResponse().SetGainNormalisationFactor(1.);
00685073 130 SwitchOffGainCorrection();
e75408ba 131 }
132 else{
133 // Assume Data
9daf22f8 134 //if(fPIDmethod == kNNpid) fPidRecal->SetGainScaleFactor(1.14);
135 //fESDpid->GetTRDResponse().SetGainNormalisationFactor(1.14);
00685073 136 SwitchOnGainCorrection();
e75408ba 137 }
e75408ba 138}
139
140//_____________________________________________________
141void AliTRDTenderSupply::ProcessEvent()
142{
143 //
144 // Reapply pid information
145 //
00685073 146 if (fTender->RunChanged()){
147 AliDebug(0, Form("AliTPCTenderSupply::ProcessEvent - Run Changed (%d)\n",fTender->GetRun()));
148 if (fGainCorrection) SetChamberGain();
51a0ce25 149 if(!fHasReferences) LoadReferences();
00685073 150 }
e75408ba 151
00685073 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);
9daf22f8 161 // Recalculate likelihoods
00685073 162 if(!(track->GetStatus() & AliESDtrack::kTRDout)) continue;
163 if(fGainCorrection) ApplyGainCorrection(track);
9daf22f8 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 }
00685073 173 }
174}
e75408ba 175
51a0ce25 176//_____________________________________________________
177void 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
00685073 221//_____________________________________________________
222void AliTRDTenderSupply::SetChamberGain(){
223 //
224 // Load Chamber Gain factors into the Tender supply
225 //
e75408ba 226
00685073 227 //find previous entry from the UserInfo
51a0ce25 228 TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree();
00685073 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
9daf22f8 254 AliCDBEntry *entry=fTender->GetCDBManager()->Get(id->GetPath(), id->GetFirstRun(), id->GetVersion());
00685073 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
9daf22f8 267 AliCDBEntry *entry=fTender->GetCDBManager()->Get(id->GetPath(), id->GetFirstRun(), id->GetVersion());
00685073 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//_____________________________________________________
300void AliTRDTenderSupply::ApplyGainCorrection(AliESDtrack * track){
e75408ba 301 //
00685073 302 // Apply new gain factors to the track
e75408ba 303 //
00685073 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//_____________________________________________________
363Bool_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");
e75408ba 387 break;
00685073 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++;
e75408ba 410 }
00685073 411 return nDet ? kTRUE : kFALSE;
e75408ba 412}
00685073 413