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 ///////////////////////////////////////////////////////////////////////////////
25 #include <TObjString.h>
29 #include <AliCDBEntry.h>
31 #include <AliCDBManager.h>
32 #include <AliTRDCalDet.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>
47 #include "AliTRDTenderSupply.h"
49 ClassImp(AliTRDTenderSupply)
51 //_____________________________________________________
52 AliTRDTenderSupply::AliTRDTenderSupply() :
56 fChamberGainOld(NULL),
57 fChamberGainNew(NULL),
58 fChamberVdriftOld(NULL),
59 fChamberVdriftNew(NULL),
63 fGainCorrection(kTRUE)
70 //_____________________________________________________
71 AliTRDTenderSupply::AliTRDTenderSupply(const char *name, const AliTender *tender) :
72 AliTenderSupply(name,tender),
75 fChamberGainOld(NULL),
76 fChamberGainNew(NULL),
77 fChamberVdriftOld(NULL),
78 fChamberVdriftNew(NULL),
82 fGainCorrection(kTRUE)
87 memset(fBadChamberID, 0, sizeof(Int_t) * kNChambers);
90 //_____________________________________________________
91 AliTRDTenderSupply::~AliTRDTenderSupply()
98 //_____________________________________________________
99 void AliTRDTenderSupply::Init()
102 // Initialise TRD tender
105 AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
107 // 1DLQ PID implemented in the AliESD object
108 fESDpid=fTender->GetESDhandler()->GetESDpid();
110 fESDpid=new AliESDpid;
111 fTender->GetESDhandler()->SetESDpid(fESDpid);
113 // Set Normalisation Factors
114 if(mgr->GetMCtruthEventHandler()){
116 SwitchOffGainCorrection();
120 SwitchOnGainCorrection();
124 //_____________________________________________________
125 void AliTRDTenderSupply::ProcessEvent()
128 // Reapply pid information
130 if (fTender->RunChanged()){
131 AliDebug(0, Form("AliTPCTenderSupply::ProcessEvent - Run Changed (%d)\n",fTender->GetRun()));
132 if (fGainCorrection) SetChamberGain();
135 fESD = fTender->GetEvent();
137 Int_t ntracks=fESD->GetNumberOfTracks();
140 // recalculate PID probabilities
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));
150 //_____________________________________________________
151 void AliTRDTenderSupply::SetChamberGain(){
153 // Load Chamber Gain factors into the Tender supply
156 //find previous entry from the UserInfo
157 TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree();
159 AliError("Tree not found in ESDhandler");
163 TList *userInfo=(TList*)tree->GetUserInfo();
165 AliError("No UserInfo found in tree");
169 TList *cdbList=(TList*)userInfo->FindObject("cdbList");
171 AliError("No cdbList found in UserInfo");
172 if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
176 TIter nextCDB(cdbList);
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());
183 AliCDBEntry *entry=fTender->GetCDBManager()->Get(*id);
185 AliError("No previous gain calibration entry found");
189 fChamberGainOld = dynamic_cast<AliTRDCalDet *>(entry->GetObject());
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());
196 AliCDBEntry *entry=fTender->GetCDBManager()->Get(*id);
198 AliError("No previous drift velocity calibration entry found");
202 fChamberVdriftOld = dynamic_cast<AliTRDCalDet *>(entry->GetObject());
204 AliDebug(1, Form("Used old Drift velocity entry: %s\n",entry->GetId().ToString().Data()));
209 // Get Latest Gain Calib Object
210 AliCDBEntry *entryNew=fTender->GetCDBManager()->Get("TRD/Calib/ChamberGainFactor",fTender->GetRun());
212 AliDebug(1, Form("Used new Gain entry: %s\n",entryNew->GetId().ToString().Data()));
213 fChamberGainNew = dynamic_cast<AliTRDCalDet *>(entryNew->GetObject());
215 AliError("No new gain calibration entry found");
217 // Also get the latest Drift Velocity calibration object
218 entryNew=fTender->GetCDBManager()->Get("TRD/Calib/ChamberVdrift",fTender->GetRun());
220 AliDebug(1, Form("Used new Drift velocity entry: %s\n",entryNew->GetId().ToString().Data()));
221 fChamberVdriftNew = dynamic_cast<AliTRDCalDet *>(entryNew->GetObject());
223 AliError("No new drift velocity calibration entry found");
225 if(!fChamberGainNew || !fChamberVdriftNew) AliError("No recent calibration found");
228 //_____________________________________________________
229 void AliTRDTenderSupply::ApplyGainCorrection(AliESDtrack * track){
231 // Apply new gain factors to the track
233 if(!fChamberGainNew || !fChamberGainOld){
234 AliError("Cannot apply gain correction.");
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
242 Bool_t applyCorrectionVdrift = kFALSE;
243 if(fChamberVdriftOld && fChamberVdriftNew) applyCorrectionVdrift = kTRUE;
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;
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]){
260 for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
261 track->SetTRDslice(0, iplane, islice);
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;
276 AliDebug(2, Form("Applying correction factor %f\n", correction));
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);
284 if(nSlices) nTrackletsPID++;
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));
291 //_____________________________________________________
292 Bool_t AliTRDTenderSupply::GetTRDchamberID(AliESDtrack * const track, Int_t *detectors) {
294 // Calculate TRD chamber ID
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;
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();
306 AliDebug(2, "No Track Params");
310 AliExternalTrackParam workparam(*trueparam); // Do calculation on working Copy
313 for(Int_t ily = 0; ily < kNPlanes; ily++){
314 if(!workparam.PropagateTo(xLayer[ily], fESD->GetMagneticField())) {
315 AliDebug(2, "Propagation failed");
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.;
323 Int_t sector = static_cast<Int_t>(trackAlpha/secAlpha);
324 Double_t etaTrack = track->Eta();
326 for(Int_t istack = 0; istack < 5; istack++){
327 if(etaTrack >= etamin[istack] && etaTrack <= etamax[istack]){
333 AliDebug(2, "Dead Area");
337 detectors[ily] = sector * kNStacks * kNPlanes + stack * kNPlanes + ily;
340 return nDet ? kTRUE : kFALSE;