1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // TRD tender: reapply pid on the fly //
21 ///////////////////////////////////////////////////////////////////////////////
24 #include <TDirectory.h>
27 #include <TObjString.h>
32 #include <AliCDBEntry.h>
34 #include <AliCDBManager.h>
35 #include <AliOADBContainer.h>
36 #include <AliTRDCalDet.h>
41 #include <AliGeomManager.h>
43 #include <AliVEvent.h>
44 #include <AliESDEvent.h>
45 #include <AliESDpid.h>
46 #include <AliESDtrack.h>
47 #include <AliESDInputHandler.h>
48 #include <AliAnalysisManager.h>
49 #include <AliTrackerBase.h>
50 #include <AliTRDPIDResponseObject.h>
51 #include <AliTRDPIDResponse.h>
52 #include "AliTRDCalChamberStatus.h"
53 #include <AliTender.h>
55 #include "AliTRDTenderSupply.h"
57 ClassImp(AliTRDTenderSupply)
59 //_____________________________________________________
60 AliTRDTenderSupply::AliTRDTenderSupply() :
64 fChamberGainOld(NULL),
65 fChamberGainNew(NULL),
66 fChamberVdriftOld(NULL),
67 fChamberVdriftNew(NULL),
68 fRunByRunCorrection(NULL),
70 fNormalizationFactor(1.),
74 fGainCorrection(kTRUE),
75 fLoadReferences(kFALSE),
76 fLoadReferencesFromCDB(kFALSE),
77 fLoadDeadChambers(kFALSE),
78 fHasReferences(kFALSE),
79 fHasNewCalibration(kTRUE),
81 fNameRunByRunCorrection()
86 memset(fSlicesForPID, 0, sizeof(UInt_t) * 2);
87 memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers);
90 //_____________________________________________________
91 AliTRDTenderSupply::AliTRDTenderSupply(const char *name, const AliTender *tender) :
92 AliTenderSupply(name,tender),
95 fChamberGainOld(NULL),
96 fChamberGainNew(NULL),
97 fChamberVdriftOld(NULL),
98 fChamberVdriftNew(NULL),
99 fRunByRunCorrection(NULL),
101 fNormalizationFactor(1.),
105 fGainCorrection(kTRUE),
106 fLoadReferences(kFALSE),
107 fLoadReferencesFromCDB(kFALSE),
108 fLoadDeadChambers(kFALSE),
109 fHasReferences(kFALSE),
110 fHasNewCalibration(kTRUE),
112 fNameRunByRunCorrection()
117 memset(fSlicesForPID, 0, sizeof(UInt_t) * 2);
118 memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers);
121 //_____________________________________________________
122 AliTRDTenderSupply::~AliTRDTenderSupply()
129 //_____________________________________________________
130 void AliTRDTenderSupply::Init()
133 // Initialise TRD tender
136 AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
138 // 1DLQ PID implemented in the AliESD object
139 fESDpid=fTender->GetESDhandler()->GetESDpid();
141 fESDpid=new AliESDpid;
142 fTender->GetESDhandler()->SetESDpid(fESDpid);
145 if(fLoadReferences && !fLoadReferencesFromCDB) LoadReferences();
146 //fESDpid->GetTRDResponse().SetGainNormalisationFactor(fNormalizationFactor);
147 fESDpid->SetTRDslicesForPID(fSlicesForPID[0], fSlicesForPID[1]);
149 if(fNameRunByRunCorrection.Length()) LoadRunByRunCorrection(fNameRunByRunCorrection.Data());
151 // Set Normalisation Factors
152 if(mgr->GetMCtruthEventHandler()){
154 //fESDpid->GetTRDResponse().SetGainNormalisationFactor(1.);
155 SwitchOffGainCorrection();
159 //if(fPIDmethod == kNNpid) fPidRecal->SetGainScaleFactor(1.14);
160 //fESDpid->GetTRDResponse().SetGainNormalisationFactor(1.14);
161 SwitchOnGainCorrection();
165 //_____________________________________________________
166 void AliTRDTenderSupply::ProcessEvent()
169 // Reapply pid information
171 if (fTender->RunChanged()){
172 AliDebug(0, Form("AliTPCTenderSupply::ProcessEvent - Run Changed (%d)\n",fTender->GetRun()));
173 if (fGainCorrection) SetChamberGain();
174 if(fLoadReferences && !fHasReferences) LoadReferences();
175 if(fLoadDeadChambers) LoadDeadChambersFromCDB();
177 if(AliGeomManager::GetGeometry()){
178 AliInfo("Geometry already loaded by other tenders");
180 if(fGeoFile) AliInfo(Form("Load geometry from file %s\n", fGeoFile));
181 else AliInfo("Load Geometry from OCDB\n");
182 AliGeomManager::LoadGeometry(fGeoFile);
187 fESD = fTender->GetEvent();
189 Int_t ntracks=fESD->GetNumberOfTracks();
192 // recalculate PID probabilities
194 Int_t detectors[kNPlanes];
195 for(Int_t itrack = 0; itrack < ntracks; itrack++){
196 for(Int_t idet = 0; idet < 5; idet++) detectors[idet] = -1;
197 AliESDtrack *track=fESD->GetTrack(itrack);
198 // Recalculate likelihoods
199 if(!(track->GetStatus() & AliESDtrack::kTRDout)) continue;
200 AliDebug(2, Form("TRD track found, gain correction: %s, Number of bad chambers: %d\n", fGainCorrection ? "Yes" : "No", fNBadChambers));
201 if(GetTRDchamberID(track, detectors)){
202 if(fGainCorrection && fHasNewCalibration) ApplyGainCorrection(track, detectors);
203 if(fNBadChambers) MaskChambers(track, detectors);
205 if(fRunByRunCorrection) ApplyRunByRunCorrection(track);
206 if(fNormalizationFactor != 1.){
207 //printf("Gain Factor: %f\n", fNormalizationFactor);
208 // Renormalize charge
209 Double_t qslice = -1;
210 for(Int_t ily = 0; ily < 6; ily++){
211 for(Int_t is = 0; is < track->GetNumberOfTRDslices(); is++){
212 qslice = track->GetTRDslice(ily, is);
213 //printf("Doing layer %d slice %d, value %f\n", ily, is, qslice);
215 qslice *= fNormalizationFactor;
216 //printf("qslice new: %f\n", qslice);
217 track->SetTRDslice(qslice, ily, is);
226 fESDpid->MakeTRDPID(fESD->GetTrack(itrack));
229 AliError("PID Method not implemented (yet)");
234 //_____________________________________________________
235 void AliTRDTenderSupply::LoadDeadChambersFromCDB(){
237 // Load Dead Chambers from the OCDB
239 AliDebug(1, "Loading Dead Chambers from the OCDB");
240 AliCDBEntry *en = fTender->GetCDBManager()->Get("TRD/Calib/ChamberStatus",fTender->GetRun());
242 AliError("Dead Chambers not in OCDB");
247 AliTRDCalChamberStatus* chamberStatus = 0;
249 chamberStatus = (AliTRDCalChamberStatus*)en->GetObject();
250 if(!chamberStatus) AliError("List with the dead chambers not found");
251 for(Int_t ichamber = 0; ichamber < 540; ichamber++) {
252 if(!chamberStatus->IsGood(ichamber)){
253 //printf("Chamber not installed %d\n",ichamber);
254 AddBadChamber(ichamber);
260 //_____________________________________________________
261 void AliTRDTenderSupply::LoadReferences(){
263 // Load Reference from the OCDB/OADB into the PID Response
265 if(fLoadReferencesFromCDB){
266 AliDebug(1, "Loading Reference Distributions from the OCDB");
267 AliCDBEntry *en = fTender->GetCDBManager()->Get("TRD/Calib/PIDLQ1D");
269 AliError("References for 1D Likelihood Method not in OCDB");
273 TObjArray *arr = dynamic_cast<TObjArray *>(en->GetObject());
274 if(!arr) AliError("List with the references not found");
276 // Get new references
279 AliTRDPIDResponseObject *ref = NULL;
281 if(!TString(o->IsA()->GetName()).CompareTo("AliTRDPIDResponseObject")){
282 ref = dynamic_cast<AliTRDPIDResponseObject *>(o);
287 fESDpid->GetTRDResponse().SetPIDResponseObject(ref);
288 fHasReferences = kTRUE;
289 AliDebug(1, "Reference distributions loaded into the PID Response");
291 AliError("References not found");
294 // Backward compatibility mode
295 AliInfo("Loading Reference Distributions from ROOT file");
296 fESDpid->GetTRDResponse().Load("$TRAIN_ROOT/util/tender/LQ1dRef_v3.root");
297 fHasReferences = kTRUE;
301 //_____________________________________________________
302 void AliTRDTenderSupply::SetChamberGain(){
304 // Load Chamber Gain factors into the Tender supply
307 //find previous entry from the UserInfo
308 TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree();
310 fHasNewCalibration = kFALSE;
311 AliError("Tree not found in ESDhandler");
315 TList *userInfo=(TList*)tree->GetUserInfo();
317 fHasNewCalibration = kFALSE;
318 AliError("No UserInfo found in tree");
322 TList *cdbList=(TList*)userInfo->FindObject("cdbList");
324 fHasNewCalibration = kFALSE;
325 AliError("No cdbList found in UserInfo");
326 if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
329 fHasNewCalibration = kTRUE;
331 TIter nextCDB(cdbList);
333 while ( (os=(TObjString*)nextCDB()) ){
334 if(os->GetString().Contains("TRD/Calib/ChamberGainFactor")){
335 // Get Old gain calibration
336 AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
338 AliCDBEntry *entry=fTender->GetCDBManager()->Get(id->GetPath(), id->GetFirstRun(), id->GetVersion());
340 AliError("No previous gain calibration entry found");
344 fChamberGainOld = dynamic_cast<AliTRDCalDet *>(entry->GetObject());
346 AliDebug(1, Form("Used old Gain entry: %s\n",entry->GetId().ToString().Data()));
347 } else if(os->GetString().Contains("TRD/Calib/ChamberVdrift")){
348 // Get Old drift velocity calibration
349 AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
351 AliCDBEntry *entry=fTender->GetCDBManager()->Get(id->GetPath(), id->GetFirstRun(), id->GetVersion());
353 AliError("No previous drift velocity calibration entry found");
357 fChamberVdriftOld = dynamic_cast<AliTRDCalDet *>(entry->GetObject());
359 AliDebug(1, Form("Used old Drift velocity entry: %s\n",entry->GetId().ToString().Data()));
364 // Get Latest Gain Calib Object
365 AliCDBEntry *entryNew=fTender->GetCDBManager()->Get("TRD/Calib/ChamberGainFactor",fTender->GetRun());
367 AliDebug(1, Form("Used new Gain entry: %s\n",entryNew->GetId().ToString().Data()));
368 fChamberGainNew = dynamic_cast<AliTRDCalDet *>(entryNew->GetObject());
370 AliError("No new gain calibration entry found");
372 // Also get the latest Drift Velocity calibration object
373 entryNew=fTender->GetCDBManager()->Get("TRD/Calib/ChamberVdrift",fTender->GetRun());
375 AliDebug(1, Form("Used new Drift velocity entry: %s\n",entryNew->GetId().ToString().Data()));
376 fChamberVdriftNew = dynamic_cast<AliTRDCalDet *>(entryNew->GetObject());
378 AliError("No new drift velocity calibration entry found");
380 if(!fChamberGainNew || !fChamberVdriftNew){
381 AliError("No recent calibration found");
382 fHasNewCalibration = kFALSE;
386 //_____________________________________________________
387 void AliTRDTenderSupply::LoadRunByRunCorrection(const char *filename){
389 // Define run by run gain correction for the charge
392 TDirectory *bkp = gDirectory;
393 TFile *in = TFile::Open(filename);
395 fRunByRunCorrection = dynamic_cast<AliOADBContainer *>(in->Get("TRDchargeCorrection"));
397 if(fRunByRunCorrection )
398 AliDebug(2, Form("OADB Container has %d runs\n", fRunByRunCorrection->GetNumberOfEntries()));
399 /* Temporarily out due to a bug in AliOADBContainer
400 fRunByRunCorrection = new AliOADBContainer("TRDchargeCorrection");
401 Int_t status = fRunByRunCorrection->InitFromFile(filename, "TRDchargeCorrection");
402 if(!status) AliDebug(1, Form("Run-dependend gain correction factors loaded from OADB file %s", filename));
404 AliDebug(1, "Failed Loading Run-dependend gain correction factors");
405 delete fRunByRunCorrection;
406 fRunByRunCorrection = NULL;
411 //_____________________________________________________
412 Bool_t AliTRDTenderSupply::IsBadChamber(Int_t chamberID){
414 // Check if the chamber id is in the list of bad chambers
416 Bool_t isBad = kFALSE;
417 for(UInt_t icam = 0; icam < fNBadChambers; icam++)
418 if(fBadChamberID[icam] == chamberID){
420 //printf("cross checking: %i \n",chamberID);
426 //_____________________________________________________
427 void AliTRDTenderSupply::ApplyGainCorrection(AliESDtrack * track, const Int_t * const chamberID){
429 // Apply new gain factors to the track
431 if(!fChamberGainNew || !fChamberGainOld){
432 AliError("Cannot apply gain correction.");
436 if(!(track->GetStatus() & AliESDtrack::kTRDout)) return;
437 Bool_t applyCorrectionVdrift = kFALSE;
438 if(fChamberVdriftOld && fChamberVdriftNew) applyCorrectionVdrift = kTRUE;
440 for(Int_t iplane = 0; iplane < kNPlanes; iplane++){
441 if(chamberID[iplane] < 0) continue;
442 if(IsBadChamber(chamberID[iplane])) continue; // Don't apply gain correction for chambers which are in the list of bad chambers
444 // Take old and new gain factor and make ratio
445 Double_t facOld = fChamberGainOld->GetValue(chamberID[iplane]);
446 Double_t facNew = fChamberGainNew->GetValue(chamberID[iplane]);
447 Double_t correction = facNew/facOld;
448 if(applyCorrectionVdrift){
449 // apply also correction for drift velocity calibration
450 Double_t vDriftOld = fChamberVdriftOld->GetValue(chamberID[iplane]);
451 Double_t vDriftNew = fChamberVdriftNew->GetValue(chamberID[iplane]);
452 correction *= vDriftNew/vDriftOld;
454 AliDebug(2, Form("Applying correction factor %f\n", correction));
455 for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
456 Double_t qslice = track->GetTRDslice(iplane, islice);
457 if(qslice <= 0.) continue;
458 track->SetTRDslice(qslice / correction, iplane, islice);
463 //_____________________________________________________
464 void AliTRDTenderSupply::ApplyRunByRunCorrection(AliESDtrack *const track) {
466 // Equalize charge distribution by applying run-by-run correction (multiplicative)
469 TVectorD *corrfactor = dynamic_cast<TVectorD *>(fRunByRunCorrection->GetObject(fTender->GetRun()));
471 AliDebug(2, "Couldn't derive gain correction factor from OADB");
474 else AliDebug(2, Form("Gain factor from OADB %f", (*corrfactor)[0]));
476 for(Int_t ily = 0; ily < kNPlanes; ily++){
477 for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
478 slice = track->GetTRDslice(ily, islice);
479 if(slice < 0.001) continue; // do not modify slices which are 0 or negative
480 slice *= (*corrfactor)[0];
481 track->SetTRDslice(slice, ily, islice);
486 //_____________________________________________________
487 void AliTRDTenderSupply::MaskChambers(AliESDtrack *const track, const Int_t * const chamberID){
489 // Mask out chambers which are in the list of bad chambers
490 // Set chamber signal to 0 and reduce the number of tracklets used for PID
492 AliDebug(2, "Masking bad chambers for TRD track");
493 Int_t nTrackletsPID = 0, nslice = 0, nTracklets = track->GetTRDntracklets();
494 Bool_t badChamber = kFALSE;
495 //Int_t nbad = 0 ; // Number of bad chambers which contain also a signal
496 //Int_t nsliceBad = 0; // Number of slices in tracklet in a bad chamber
497 for(Int_t iplane = 0; iplane < kNPlanes; iplane++){
499 nslice = 0; //nsliceBad = 0;
500 if(IsBadChamber(chamberID[iplane])) badChamber = kTRUE;
501 for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
503 //if(track->GetTRDslice(iplane, islice)) nsliceBad++;
504 track->SetTRDslice(-1, iplane, islice);
505 } else if(track->GetTRDslice(iplane, islice) > 0.001) nslice++;
507 //if(nsliceBad) nbad++;
508 if(nslice > 0) nTrackletsPID++;
510 //if(nbad) track->SetTRDncls(track->GetTRDncls() - 20 * nbad); // subtract mean number of clusters per tracklet for bad tracklets
511 // Use nTrackletsPID to indicate the number of tracklets from good
512 // chambers so they are used for the PID
513 track->SetTRDntracklets(nTrackletsPID | (nTracklets << 3));
516 //_____________________________________________________
517 Bool_t AliTRDTenderSupply::GetTRDchamberID(AliESDtrack * const track, Int_t *detectors) {
519 // Calculate TRD chamber ID
521 Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
522 if(p < fPthreshold) return kFALSE; // Apply low momentum cutoff
524 Double_t xLayer[kNPlanes] = {300.2, 312.8, 325.4, 338., 350.6, 363.2};
525 Double_t etamin[kNStacks] = {0.536, 0.157, -0.145, -0.527,-0.851};
526 Double_t etamax[kNStacks] = {0.851, 0.527, 0.145, -0.157,-0.536};
527 for(Int_t ily = 0; ily < kNPlanes; ily++) detectors[ily] = -1;
529 const AliExternalTrackParam *trueparam = NULL;
530 if(track->GetOuterParam()) trueparam = track->GetOuterParam();
531 else if(track->GetTPCInnerParam()) trueparam = track->GetTPCInnerParam();
532 else if(track->GetInnerParam()) trueparam = track->GetInnerParam();
534 AliDebug(2, "No Track Params");
538 AliExternalTrackParam workparam(*trueparam); // Do calculation on working Copy
541 for(Int_t ily = 0; ily < kNPlanes; ily++){
542 if(!AliTrackerBase::PropagateTrackToBxByBz(&workparam, xLayer[ily], 0.139, 100)){ // Assuming the pion mass
543 AliDebug(2, "Propagation failed");
546 workparam.GetXYZ(pos);
547 Double_t trackAlpha = TMath::ATan2(pos[1], pos[0]);
549 // Compare to simple propagation without magnetic field
550 AliExternalTrackParam workparam1(*trueparam); // Do calculation on working Copy
552 if(!workparam1.PropagateTo(xLayer[ily], fESD->GetMagneticField())) {
553 AliDebug(2, "Propagation failed");
556 workparam1.GetXYZ(pos1);
557 Double_t trackAlpha1 = TMath::ATan2(pos1[1], pos1[0]);
558 AliDebug(2, Form("Alpha: Old %f, New %f, diff %f", trackAlpha1, trackAlpha, trackAlpha-trackAlpha1));
560 if(trackAlpha < 0) trackAlpha = 2 * TMath::Pi() + trackAlpha;
561 Double_t secAlpha = 2 * TMath::Pi() / 18.;
563 Int_t sector = static_cast<Int_t>(trackAlpha/secAlpha);
564 Double_t etaTrack = track->Eta();
566 for(Int_t istack = 0; istack < 5; istack++){
567 if(etaTrack >= etamin[istack] && etaTrack <= etamax[istack]){
573 AliDebug(2, "Dead Area");
577 detectors[ily] = sector * kNStacks * kNPlanes + stack * kNPlanes + ily;
580 return nDet ? kTRUE : kFALSE;