9 #include "AliMCParticle.h"
10 //#include "AliStack.h"
12 #include "AliAnalysisTask.h"
13 #include "AliAnalysisManager.h"
15 #include "AliESDEvent.h"
16 #include "AliMCEvent.h"
17 #include "AliESDInputHandler.h"
18 #include "AliInputEventHandler.h"
20 #include "AliVVertex.h"
21 #include "AliAnalysisFilter.h"
23 #include "AliPIDResponse.h"
24 #include "AliTPCPIDResponse.h"
25 //#include "AliTPCParamSR.h"
27 #include "AliTPCPIDEtaTree.h"
30 This task determines the eta dependence of the TPC signal.
31 For this purpose, only tracks fulfilling some standard quality cuts are taken into account.
32 The obtained data can be used to derive the functional behaviour of the eta dependence.
33 Such a function can be plugged into this task to correct for the eta dependence and to see
34 if there is then really no eta dependence left.
36 Class written by Benjamin Hess.
37 Contact: bhess@cern.ch
40 ClassImp(AliTPCPIDEtaTree)
42 //________________________________________________________________________
43 AliTPCPIDEtaTree::AliTPCPIDEtaTree()
45 , fStoreMultiplicity(kFALSE)
46 , fStoreNumOfSubthresholdclusters(kFALSE)
47 , fStoreNumClustersInActiveVolume(kFALSE)
48 , fDoAdditionalQA(kFALSE)
59 , fTPCsignalNsubthreshold(0)
60 , fNumTPCClustersInActiveVolume(0)
63 , fCorrectdEdxEtaDependence(kFALSE)
64 , fCorrectdEdxMultiplicityDependence(kFALSE)
67 , fOutputContainer(0x0)
69 , fhMultiplicityQA(0x0)
71 // default Constructor
74 //________________________________________________________________________
75 AliTPCPIDEtaTree::AliTPCPIDEtaTree(const char *name)
77 , fStoreMultiplicity(kFALSE)
78 , fStoreNumOfSubthresholdclusters(kFALSE)
79 , fStoreNumClustersInActiveVolume(kFALSE)
80 , fDoAdditionalQA(kFALSE)
91 , fTPCsignalNsubthreshold(0)
92 , fNumTPCClustersInActiveVolume(0)
95 , fCorrectdEdxEtaDependence(kFALSE)
96 , fCorrectdEdxMultiplicityDependence(kFALSE)
99 , fOutputContainer(0x0)
101 , fhMultiplicityQA(0x0)
105 // Define input and output slots here
106 DefineInput(0, TChain::Class());
107 DefineOutput(1, TTree::Class());
108 DefineOutput(2, TTree::Class());
109 DefineOutput(3, TObjArray::Class());
113 //________________________________________________________________________
114 AliTPCPIDEtaTree::~AliTPCPIDEtaTree()
121 delete fOutputContainer;
122 fOutputContainer = 0x0;
126 //________________________________________________________________________
127 void AliTPCPIDEtaTree::UserCreateOutputObjects()
132 AliTPCPIDBase::UserCreateOutputObjects();
134 if (fDoAdditionalQA) {
137 fOutputContainer = new TObjArray(2);
138 fOutputContainer->SetName(GetName());
139 fOutputContainer->SetOwner(kTRUE);
141 const Int_t nBins = 6;
142 // p_vert, p_TPC, eta, nSigmaTOF, beta, multiplicity
143 Int_t bins[nBins] = { 48, 48, 9, 30, 110, 4 };
144 Double_t xmin[nBins] = { 0.3, 0.3, -0.9, -5.0, 0.0, 0. };
145 Double_t xmax[nBins] = { 2.7, 2.7, 0.9, 5.0, 1.1, 3200 };
148 fhTOFqa = new THnSparseI("hTOFqa","", nBins, bins, xmin, xmax);
149 fhTOFqa->GetAxis(0)->SetTitle("p_{Vertex} (GeV/c)");
150 fhTOFqa->GetAxis(1)->SetTitle("p_{TPC_inner} (GeV/c)");
151 fhTOFqa->GetAxis(2)->SetTitle("#eta");
152 fhTOFqa->GetAxis(3)->SetTitle("n #sigma_{p} TOF");
153 fhTOFqa->GetAxis(4)->SetTitle("#beta");
154 fhTOFqa->GetAxis(5)->SetTitle("Multiplicity");
156 fOutputContainer->Add(fhTOFqa);
160 fhMultiplicityQA = new TH2I("hMultiplicityQA", "Multiplicity check; Contributors to primary vertex per event; Total number of tracks per Event",
161 120, 0, 6000, 400, 0, 20000);
162 fOutputContainer->Add(fhMultiplicityQA);
165 fOutputContainer = new TObjArray(1);
166 fOutputContainer->SetName(GetName());
167 fOutputContainer->SetOwner(kTRUE);
172 fTree = new TTree("fTree", "Tree for analysis of #eta dependence of TPC signal");
173 fTree->Branch("pTPC", &fPtpc);
174 //fTree->Branch("pT", &fPt);
175 fTree->Branch("dEdx", &fDeDx);
176 fTree->Branch("dEdxExpected", &fDeDxExpected);
177 fTree->Branch("tanTheta", &fTanTheta);
178 //fTree->Branch("sinAlpha", &fSinAlpha);
179 //fTree->Branch("y", &fY);
180 //TODO fTree->Branch("phiPrime", &fPhiPrime);
181 fTree->Branch("tpcSignalN", &fTPCsignalN);
183 if (fStoreNumOfSubthresholdclusters)
184 fTree->Branch("tpcSignalNsubthreshold", &fTPCsignalNsubthreshold);
186 if (fStoreNumClustersInActiveVolume)
187 fTree->Branch("numTPCClustersInActiveVolume", &fNumTPCClustersInActiveVolume);
189 fTree->Branch("pidType", &fPIDtype);
191 if (fStoreMultiplicity) {
192 fTree->Branch("fMultiplicity", &fMultiplicity);
197 fTreePions = new TTree("fTreePions", "Tree for analysis of #eta dependence of TPC signal - Pions");
198 fTreePions->Branch("pTPC", &fPtpc);
199 fTreePions->Branch("pT", &fPt);
200 fTreePions->Branch("dEdx", &fDeDx);
201 fTreePions->Branch("dEdxExpected", &fDeDxExpected);
202 fTreePions->Branch("tanTheta", &fTanTheta);
203 fTreePions->Branch("tpcSignalN", &fTPCsignalN);
205 if (fStoreNumOfSubthresholdclusters)
206 fTreePions->Branch("tpcSignalNsubthreshold", &fTPCsignalNsubthreshold);
208 if (fStoreMultiplicity) {
209 fTreePions->Branch("fMultiplicity", &fMultiplicity);
213 PostData(2, fTreePions);
214 PostData(3, fOutputContainer);
218 //________________________________________________________________________
219 void AliTPCPIDEtaTree::UserExec(Option_t *)
222 // Called for each event
224 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
226 Printf("ERROR: fESD not available");
230 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
232 if (!fPIDResponse || !fV0KineCuts)
235 if (!GetVertexIsOk(fESD))
238 const AliVVertex* primaryVertex = fESD->GetPrimaryVertexTracks();
241 //fMultiplicity = primaryVertex->GetNContributors();
244 if(primaryVertex->GetNContributors() <= 0)
247 fMultiplicity = fESD->GetNumberOfTracks();
249 if (fDoAdditionalQA) {
250 Int_t nTotTracks = fESD->GetNumberOfTracks();
251 fhMultiplicityQA->Fill(primaryVertex->GetNContributors(), nTotTracks);
254 Double_t magField = fESD->GetMagneticField();
256 // Fill V0 arrays with V0s
262 // Track loop to fill a Train spectrum
263 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
264 Bool_t isPr = kFALSE;
265 Bool_t isPi = kFALSE;
267 AliESDtrack* track = fESD->GetTrack(iTracks);
269 Printf("ERROR: Could not receive track %d", iTracks);
273 if (TMath::Abs(track->Eta()) > fEtaCut)
276 // Do not distinguish positively and negatively charged V0's
277 Char_t v0tagAllCharges = TMath::Abs(GetV0tag(iTracks));
278 if (v0tagAllCharges == -99) {
279 AliError(Form("Problem with V0 tag list (requested V0 track for track %d from %d (list status %d))!", iTracks, fESD->GetNumberOfTracks(),
284 Bool_t isV0prNotMC = (v0tagAllCharges == 16) && !fMC; // Only accept V0 protons for data, not for MC
285 Bool_t isV0piNotMC = (v0tagAllCharges == 15) && !fMC; // Only accept V0 pions for data, not for MC
287 Bool_t isV0NotMC = isV0prNotMC || isV0piNotMC;
290 if(!isV0NotMC && fTrackFilter && !fTrackFilter->IsSelected(track))
293 // Note: For V0's, the cut on ncl can be done via the tree (value is stored). One should not cut on geometry
294 // for V0's since they have different topology
295 if (!isV0NotMC && GetUseTPCCutMIGeo()) {
296 if (!TPCCutMIGeo(track, InputEvent()))
301 fPtpc = track->GetTPCmomentum();
307 fPhiPrime = GetPhiPrime(track->Phi(), magField, track->Charge());
309 Int_t label = track->GetLabel();
311 AliMCParticle* mcTrack = 0x0;
315 continue; // If MC is available, reject tracks with negative ESD label
316 mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
318 Printf("ERROR: Could not receive mcTrack with label %d for track %d", label, iTracks);
322 /*// Only accept MC primaries
323 if (!fMC->Stack()->IsPhysicalPrimary(TMath::Abs(label))) {
329 Int_t pdg = mcTrack->PdgCode();
331 if (TMath::Abs(pdg) == 2212) // Proton
333 else if ((pdg == 111 || TMath::Abs(pdg) == 211) && fPtpc <= fPtpcPionCut) // Pion below desired momentum threshold
336 continue; // Only take protons and pions
340 if (pdg == 111 || TMath::Abs(pdg) == 211) {//Pion
343 else if (TMath::Abs(pdg) == 311 || TMath::Abs(pdg) == 321) {//Kaon
346 else if (TMath::Abs(pdg) == 2212) {//Proton
349 else if (TMath::Abs(pdg) == 11) {//Electron
352 else if (TMath::Abs(pdg) == 13) {//Muon
357 fDeDx = track->GetTPCsignal();
359 UInt_t status = track->GetStatus();
360 Bool_t hasTOFout = status&AliESDtrack::kTOFout;
361 Bool_t hasTOFtime = status&AliESDtrack::kTIME;
362 Bool_t hasTOF = kFALSE;
363 if (hasTOFout && hasTOFtime)
365 Float_t length = track->GetIntegratedLength();
366 // Length check only for primaries, not for V0's!
367 if (length < 350 && !isV0NotMC)
371 // Note: Normally, the track cuts include a cut on primaries. Therefore, if the particle is a V0, it would usually
372 // NOT be found via the primary cuts. This means that the PIDtype describes more or less disjoint sets
378 fPIDtype = kV0idNoTOF;
381 if (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) < 3.0)
382 fPIDtype = kV0idPlusTOFaccepted;
384 fPIDtype = kV0idPlusTOFrejected;
387 else if (isV0piNotMC) {
389 if (fPtpc > fPtpcPionCut || fDeDx > 140.) // Reject pions above desired momentum threshold and also reject too high dEdx
395 fPIDtype = kV0idNoTOF;
398 if (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion)) < 3.0)
399 fPIDtype = kV0idPlusTOFaccepted;
401 fPIDtype = kV0idPlusTOFrejected;
406 isPr = kTRUE; // If particle is accepted, it is a proton (for pions, only V0s are used)
408 if (fPtpc < 4.0 && //TODO was 2.7 // Do not accept non-V0s above this threshold -> High contamination!!!
409 (fDeDx >= 50. / TMath::Power(fPtpc, 1.3))) {// Pattern recognition instead of TPC cut to be ~independent of old TPC expected dEdx
414 else if (hasTOF && TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) < 3.0) {
415 fPIDtype = kTPCandTOFid;
418 continue; // Reject particle
421 continue; // Reject particle
425 if (fDoAdditionalQA) {
426 Double_t tofTime = track->GetTOFsignal();//in ps
427 Double_t tof = tofTime * 1E-3; // ns, average T0 fill subtracted, no info from T0detector
428 Double_t length2 = track->GetIntegratedLength();
429 Double_t c = TMath::C() * 1.E-9;// m/ns
430 length2 = length2 * 0.01; // in meters
432 Double_t beta = length2 / tof;
434 Double_t nSigmaTOF = hasTOF ? fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton) : 999;
436 // p_vert, p_TPC, eta, nSigmaTOF, beta, multiplicity
437 Double_t entry[6] = { track->P(), fPtpc, track->Eta(), nSigmaTOF, beta, fMultiplicity };
438 fhTOFqa->Fill(entry);
441 // Prepare entry for tree (some quantities have already been set)
442 // Turn eta correction off -> Important here is the pure spline value and the selection via cuts. The correction
443 // can be re-done manually, if needed. But it cannot be undone!
446 fDeDxExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
447 fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
449 fDeDxExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
450 fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
453 fTanTheta = track->GetInnerParam()->GetTgl();
454 //fSinAlpha = track->GetInnerParam()->GetSnp();
455 //fY = track->GetInnerParam()->GetY();
456 fTPCsignalN = track->GetTPCsignalN();
458 if (fStoreNumOfSubthresholdclusters) {
459 const AliTPCdEdxInfo *clsInfo = track->GetTPCdEdxInfo();
462 Double32_t signal[4] = {0}; // 0: IROC only; 1: OROC short padas; 2: OROC long pads; 3:OROC combined
463 Char_t ncl[3] = {0}; // clusters above threshold; 0: IROC; 1: OROC short; 2: OROC long
464 Char_t nrows[3] = {0}; // clusters above and below threshold; 0: IROC; 1: OROC short; 2: OROC long
466 clsInfo->GetTPCSignalRegionInfo(signal, ncl, nrows);
468 fTPCsignalNsubthreshold = nrows[0] + nrows[1] + nrows[2] - ncl[0] - ncl[1] - ncl[2];
471 fTPCsignalNsubthreshold = 200;// Set to invalid value
475 if (fStoreNumClustersInActiveVolume)
476 fNumTPCClustersInActiveVolume = track->GetLengthInActiveZone(1, 1.8, 220, magField);
478 fNumTPCClustersInActiveVolume = 200.;//Set to invalid value
483 const Double_t tanTheta = track->GetInnerParam()->GetTgl();
485 // Constant in formula for B in kGauss (factor 0.1 to convert B from Tesla to kGauss),
486 // pT in GeV/c (factor c*1e-9 to arrive at GeV/c) and curvature in 1/cm (factor 0.01 to get from m to cm)
487 const Double_t constant = TMath::C()* 1e-9 * 0.1 * 0.01;
488 const Double_t curvature = magField * constant / track->Pt(); // in 1./cm
490 //const Double_t angleIROC = TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(0, par.GetNRow(0) / 2.) * curvature) * 0.5, 1.));
491 //const Double_t angleOROC = TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(36, par.GetNRow(36) / 2.) * curvature) * 0.5, 1.));
493 Double_t averageddzdr = 0.;
496 for (Double_t r = 85; r < 245; r++) {
497 Double_t sinPhiLocal = TMath::Abs(r*curvature*0.5);
499 // Cut on |sin(phi)| as during reco
500 if (TMath::Abs(sinPhiLocal) <= 0.95) {
501 const Double_t phiLocal = TMath::ASin(sinPhiLocal);
502 const Double_t tanPhiLocal = TMath::Tan(phiLocal);
504 averageddzdr += tanTheta * TMath::Sqrt(1. + tanPhiLocal * tanPhiLocal);
510 averageddzdr /= nParts;
512 AliError("Problems during determination of dz/dr. Skipping track!");
517 //printf("padRow 0 / 36: %f / %f\n", par.GetPadRowRadii(0, par.GetNRow(0) / 2.), par.GetPadRowRadii(36, par.GetNRow(36) / 2.));
518 //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\nIROC/OROC/averageOld/average: %f / %f / %f / //%f\ntanThetaGlobalFromTheta/tanTheta/dzdr: %f / %f / %f\n\n",
519 // track->Pt(), constant, magField, 1./curvature, angleIROC, angleOROC, angleAverage, averageAngle, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, dzdr);
520 //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\ntanThetaGlobalFromTheta/tanTheta/Averageddzdr: %f / %f / %f\n\n",
521 // track->Pt(), constant, magField, 1./curvature, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, averageddzdr);
524 fTanTheta = averageddzdr;
535 PostData(2, fTreePions);
538 PostData(3, fOutputContainer);
540 // Clear the V0 PID arrays
544 //________________________________________________________________________
545 void AliTPCPIDEtaTree::Terminate(const Option_t *)
547 // Called once at the end of the query