1 //---------------------------------------------------------------------------------
2 // The example of usage of AliKFParticle & AliKFVertex classes for V0 analysis
4 // @author S.Gorbunov, I.Kisel
8 // The AliKFParticleTest macro contains a toy V0 finder for ESD tracks.
9 // At the first step, event primary vertex is reconstructed.
10 // At the second step, ideal PID hypothesis are assigned to all the particles.
11 // At the third step, V0 candidates are constructed for each pair
12 // of positive and negative particles.
13 // V0 candidate considered as good when it passes Chi^2 cut and
14 // it is placed >= 3 Sigma away from the primary vertex.
15 // Invariant mass distribution for all good V0 candidates is plotted.
17 // -= Copyright © ALICE HLT Group =-
18 //_________________________________________________________________________________
21 #ifndef ALIKFPARTICLETEST_H
22 #define ALIKFPARTICLETEST_H
27 #include "AliAnalysisTaskRL.h"
29 class AliKFParticleTest: public AliAnalysisTaskRL {
31 AliKFParticleTest(const char *name);
32 virtual ~AliKFParticleTest() {}
34 virtual void ConnectInputData(Option_t *);
35 virtual void CreateOutputObjects();
36 virtual void Exec(Option_t *option);
37 virtual void Terminate(Option_t *);
40 AliESD *fESD; //ESD object
42 TH1D *fHistoMassSignal;
46 TObjArray * fOutputContainer; // ! output data container
48 ClassDef(AliKFParticleTest, 0); // example of analysis
52 #include "Riostream.h"
56 #include "TParticle.h"
61 #include "AliTracker.h"
63 #include "TDatabasePDG.h"
64 #include "AliKFParticle.h"
65 #include "AliKFVertex.h"
68 ClassImp(AliKFParticleTest);
71 //________________________________________________________________________
72 AliKFParticleTest::AliKFParticleTest(const char *name) :AliAnalysisTaskRL(name,""), fESD(0) {
74 // Input slot #0 works with an Ntuple
75 DefineInput(0, TChain::Class());
76 // Output slot #0 writes into a TH1 container
77 DefineOutput(0, TObjArray::Class());
80 //___________________________________________________________________________
81 void AliKFParticleTest::ConnectInputData(Option_t *) {
82 // Initialize branches.
83 printf(" ConnectInputData of task %s\n", GetName());
85 char ** address = (char **)GetBranchAddress(0, "ESD");
86 if (address) fESD = (AliESD*)(*address);
89 SetBranchAddress(0, "ESD", &fESD);
95 //___________________________________________________________________________
96 void AliKFParticleTest::CreateOutputObjects() {
97 printf(" CreateOutputObjects of task %s\n", GetName());
99 if ( !gROOT->IsBatch() ) fCanvas = new TCanvas();
102 fHistoMassAll = new TH1D("massAll","AliKFParticle test", 500,0,3);
103 fHistoMassAll->SetXTitle("V^{0} invariant mass [GeV]");
104 fHistoMassAll->SetLineColor(kGreen);
105 fHistoMassAll->SetFillColor(kGreen);
107 fHistoMassSignal = new TH1D("massSignal","V^{0} signal", 500,0,3);
108 fHistoMassSignal->SetXTitle("V^{0} invariant mass [GeV]");
109 fHistoMassSignal->SetLineColor(kBlue);
110 fHistoMassSignal->SetFillColor(kBlue);
112 fOutputContainer = new TObjArray(1);
113 fOutputContainer->SetName(GetName()) ;
114 fOutputContainer->Add(fHistoMassAll);
115 fOutputContainer->Add(fHistoMassSignal);
119 void AliKFParticleTest::DrawV0()
121 //* Draw the invariant mass histogram
123 if ( gROOT->IsBatch() ) return;
125 const Int_t histoPDG [4]= {22,310,3122,421};
126 const Char_t* histoName[4]= {"#gamma","K^{0}_{s}","#Lambda","D^{0}"};
129 if( !fCanvas ) return;
133 for( Int_t i=0; i<4; i++ ){
134 histoText[i].SetTextColor(kBlue);
137 fHistoMassAll->Draw();
138 fHistoMassSignal->Draw("same");
140 for( Int_t i=0; i<4; i++ ){
141 Double_t mass = TDatabasePDG::Instance()->GetParticle(histoPDG[i])->Mass();
142 Int_t bin = fHistoMassSignal->FindBin(mass) -5;
145 for( Int_t j=bin; j<bin+10; j++ )
146 if( max<fHistoMassSignal->GetBinContent(j) ) max = fHistoMassSignal->GetBinContent(j);
147 if(max>0) histoText[i].DrawLatex( mass+.05, max, histoName[i] );
152 //________________________________________________________________________
153 void AliKFParticleTest::Exec(Option_t *) {
155 static Int_t iEvent = 0;
156 if( ++iEvent%100 ==0 ) cout<<"Event "<<iEvent<<endl;
159 TTree *tinput = (TTree*)GetInputData(0);
160 Long64_t ientry = tinput->GetReadEntry();
161 if (AliAnalysisTaskRL::GetEntry(ientry) == kFALSE) {
162 printf("Couldn't get event from the runLoader\n");
167 Double_t Bz = fESD->GetMagneticField();
168 AliKFParticle::SetField( Bz );
170 if (ientry==0) Notify ();
172 //cout <<"Event "<<ientry<<endl;
174 AliStack* stack = GetStack();
176 AliDebug(AliLog::kError, "Stack not available");
185 AliKFParticle fParticle; //* assigned KFParticle
186 Bool_t fPrimUsedFlag; //* flag shows that the particle was used for primary vertex fit
187 Bool_t fOK; //* is the track good enough
188 Int_t mcPDG; //* Monte Carlo PDG code
189 Int_t mcMotherID; //* Monte Carlo ID of its mother
192 Int_t nESDTracks=fESD->GetNumberOfTracks();
193 if( nESDTracks>1000 ) nESDTracks=1000;
195 TESDTrackInfo ESDTrackInfo[1000]; //* parallel to ESD tracks
197 //* Fill ESDTrackInfo array
199 for (Int_t iTr=0; iTr<nESDTracks; iTr++)
201 TESDTrackInfo &info = ESDTrackInfo[iTr];
203 info.fPrimUsedFlag = 0;
205 info.mcMotherID = -1;
207 //* track quality check
209 AliESDtrack *pTrack = fESD->GetTrack(iTr);
210 if( !pTrack ) continue;
211 if (pTrack->GetKinkIndex(0)>0) continue;
212 if ( !( pTrack->GetStatus()&AliESDtrack::kITSrefit ) ) continue;
214 if( pTrack->GetITSclusters(indi) <5 ) continue;
215 Int_t PDG = ( pTrack->Get1Pt() <0 ) ?321 :211;
219 Int_t mcID = TMath::Abs(pTrack->GetLabel());
220 TParticle * part = stack->Particle(TMath::Abs(mcID));
222 info.mcPDG = part->GetPdgCode();
224 if( mcID>=0 ) info.mcMotherID = part->GetFirstMother();
228 //* Construct KFParticle for the track
230 info.fParticle = AliKFParticle( *pTrack, PDG );
234 //* Find event primary vertex
238 const AliKFParticle * vSelected[1000]; //* Selected particles for vertex fit
239 Int_t vIndex[1000]; //* Indices of selected particles
240 Bool_t vFlag[1000]; //* Flags returned by the vertex finder
243 for( Int_t i = 0; i<nESDTracks; i++){
244 if(ESDTrackInfo[i].fOK ){
245 vSelected[nSelected] = &(ESDTrackInfo[i].fParticle);
246 vIndex[nSelected] = i;
250 primVtx.ConstructPrimaryVertex( vSelected, nSelected, vFlag, 3. );
251 for( Int_t i = 0; i<nSelected; i++){
252 if( vFlag[i] ) ESDTrackInfo[vIndex[i]].fPrimUsedFlag = 1;
254 if( primVtx.GetNDF() <1 ) return; //* Less then two tracks in primary vertex
259 for( Int_t iTr = 0; iTr<nESDTracks; iTr++ ){ //* first daughter
261 TESDTrackInfo &info = ESDTrackInfo[iTr];
262 if( !info.fOK ) continue;
264 for( Int_t jTr = iTr+1; jTr<nESDTracks; jTr++ ){ //* second daughter
265 TESDTrackInfo &jnfo = ESDTrackInfo[jTr];
266 if( !jnfo.fOK ) continue;
268 //* check for different charge
270 if( info.fParticle.GetQ() == jnfo.fParticle.GetQ() ) continue;
272 //* construct V0 mother
274 AliKFParticle V0( info.fParticle, jnfo.fParticle );
278 if( V0.GetNDF()<1 ) continue;
279 if( TMath::Sqrt(TMath::Abs(V0.GetChi2()/V0.GetNDF())) >3. ) continue;
281 //* subtruct daughters from primary vertex
283 AliKFVertex primVtxCopy = primVtx;
285 if( info.fPrimUsedFlag ) primVtxCopy -= info.fParticle;
286 if( jnfo.fPrimUsedFlag ) primVtxCopy -= jnfo.fParticle;
288 //* Check V0 Chi^2 deviation from primary vertex
290 if( V0.GetDeviationFromVertex( primVtxCopy ) >3. ) continue;
292 //* Add V0 to primary vertex to improve the primary vertex resolution
296 //* Set production vertex for V0
298 V0.SetProductionVertex( primVtxCopy );
300 //* Check chi^2 for a case
302 if( TMath::Sqrt( TMath::Abs(V0.GetChi2()/V0.GetNDF()) >3. )) continue;
304 //* Get V0 decay length with estimated error
306 Double_t length, sigmaLength;
307 if( V0.GetDecayLength( length, sigmaLength ) ) continue;
309 //* Reject V0 if it decays too close to the primary vertex
311 if( length <3.*sigmaLength ) continue;
313 //* Get V0 invariant mass and plot it
315 Double_t mass, sigmaMass;
316 if( V0.GetMass( mass, sigmaMass ) ) continue;
317 fHistoMassAll->Fill(mass);
319 //* Fill signal histograms using MC information
321 if( info.mcMotherID==jnfo.mcMotherID && info.mcMotherID>=0 ){
322 TParticle *mother = stack->Particle(info.mcMotherID);
323 if( mother && TMath::Abs(mother->GetPdgCode())!=21 ){
324 if( mother->GetNDaughters()==2 ){
325 fHistoMassSignal->Fill(mass);
327 cout<<"PDG V0,pi,pj, ndaughters, mc mass, reco mass = "<<mother->GetPdgCode()<<","<<info.mcPDG<<","<<jnfo.mcPDG<<", "
328 << mother->GetNDaughters()<<", "<<mother->GetMass()<<", "<<mass<<endl;
333 if( iEvent %1000 == 0 || (iEvent ==200)) DrawV0();
335 // Post final data. It will be written to a file with option "RECREATE"
336 PostData(0, fOutputContainer);
340 //________________________________________________________________________
341 void AliKFParticleTest::Terminate(Option_t *)
344 fOutputContainer = (TObjArray*)GetOutputData(0);
345 //fHistoMassAll=(TH1D*)fOutputContainer->At(0) ;