752061a83e67f6975e911a6628589cc3a2f39350
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDcheckESD.cxx
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 #include <TClonesArray.h>
18 #include <TObjArray.h>
19 #include <TObject.h>
20 #include <TH2I.h>
21 #include <TFile.h>
22 #include <TTree.h>
23 #include <TROOT.h>
24 #include <TChain.h>
25 #include <TParticle.h>
26
27 #include "AliLog.h"
28 #include "AliAnalysisManager.h"
29 #include "AliESDEvent.h"
30 #include "AliMCEvent.h"
31 #include "AliESDInputHandler.h"
32 #include "AliMCEventHandler.h"
33
34 #include "AliESDtrack.h"
35 #include "AliMCParticle.h"
36 #include "AliPID.h"
37 #include "AliStack.h"
38 #include "AliTrackReference.h"
39 #include "AliTRDgeometry.h"
40
41 #include "AliTRDcheckESD.h"
42
43 ClassImp(AliTRDcheckESD)
44
45 const Float_t AliTRDcheckESD::xTPC = 290.;
46 const Float_t AliTRDcheckESD::xTOF = 365.;
47
48 //____________________________________________________________________
49 AliTRDcheckESD::AliTRDcheckESD():
50   AliAnalysisTask("ESDchecker", "TRD checker @ ESD level")
51   ,fStatus(0)
52   ,fESD(0x0)
53   ,fMC(0x0)
54   ,fHistos(0x0)
55 {
56   //
57   // Default constructor
58   //
59
60   DefineInput(0, TChain::Class());
61   DefineOutput(0, TObjArray::Class());
62 }
63
64 //____________________________________________________________________
65 AliTRDcheckESD::~AliTRDcheckESD()
66 {
67   if(fHistos){
68     //fHistos->Delete();
69     delete fHistos;
70   }
71 }
72
73 //____________________________________________________________________
74 void AliTRDcheckESD::ConnectInputData(Option_t *)
75 {
76   //
77   // Link the Input Data
78   //
79   TTree *tree = dynamic_cast<TChain*>(GetInputData(0));
80   if(tree) tree->SetBranchStatus("Tracks", 1);
81
82   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
83   fESD = esdH ? esdH->GetEvent() : 0x0;
84
85   if(!HasMC()) return;
86   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
87   fMC = mcH ? mcH->MCEvent() : 0x0;
88 }
89
90 //____________________________________________________________________
91 void AliTRDcheckESD::CreateOutputObjects()
92 {       
93   //
94   // Create Output Containers (TObjectArray containing 1D histograms)
95   //
96   OpenFile(0, "RECREATE");  
97   fHistos = new TObjArray(5);
98   //fHistos->SetOwner(kTRUE);
99   
100   TH1 *h = 0x0;
101
102   // clusters per tracklet
103   if(!(h = (TH1I*)gROOT->FindObject("hNCl"))){
104     h = new TH1I("hNCl", "Clusters per TRD track", 100, 0., 200.);
105     h->GetXaxis()->SetTitle("N_{cl}^{TRD}");
106     h->GetYaxis()->SetTitle("entries");
107   } else h->Reset();
108   fHistos->AddAt(h, kNCl);
109
110   // TPC out
111   if(!(h = (TH2I*)gROOT->FindObject("hTRDstat"))){
112     h = new TH2I("hTRDstat", "TRD status bits", 100, 0., 20., 4, -.5, 3.5);
113     h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
114     h->GetYaxis()->SetTitle("status bits");
115     h->GetZaxis()->SetTitle("entries");
116   } else h->Reset();
117   fHistos->AddAt(h, kTRDstat);
118 }
119
120 //____________________________________________________________________
121 void AliTRDcheckESD::Exec(Option_t *){
122   //
123   // Run the Analysis
124   //
125   if(!fESD){
126     AliError("ESD not found");
127     return;
128   }
129
130   // Get MC information if available
131   AliStack * fStack = 0x0;
132   if(HasMC() && !fMC){ 
133     AliWarning("Monte Carlo Event not available");
134     SetMC(kFALSE);
135   } else {
136     if(!(fStack = fMC->Stack())){
137       AliWarning("Cannot get the Monte Carlo Stack");
138       SetMC(kFALSE);
139     }
140   }
141   Bool_t TPCout(0), TRDin(0), TRDout(0), TRDpid(0);
142
143
144   Int_t nTRD = 0, nTPC = 0;
145   //Int_t nTracks = fESD->GetNumberOfTracks();
146   AliESDtrack *esdTrack = 0x0;
147   for(Int_t itrk = 0; itrk < fESD->GetNumberOfTracks(); itrk++){
148     TPCout=0;TRDin=0;TRDout=0;TRDpid=0;
149     esdTrack = fESD->GetTrack(itrk);
150     if(esdTrack->GetNcls(1)) nTPC++;
151     if(esdTrack->GetNcls(2)) nTRD++;
152
153     // track status
154     ULong_t status = esdTrack->GetStatus();
155
156     // TRD PID
157     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
158     // pid quality
159     esdTrack->GetTRDpidQuality();
160     // kink index
161     esdTrack->GetKinkIndex(0);
162     // TPC clusters
163     esdTrack->GetNcls(1);
164
165     // look at external track param
166     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
167     Double_t xyz[3];
168     if(op){
169       op->GetXYZ(xyz);
170       op->Global2LocalPosition(xyz, op->GetAlpha());
171       //printf("op @ X[%7.3f]\n", xyz[0]);
172     }
173
174     // read MC info
175     if(!HasMC()) continue;
176
177     Int_t fLabel = esdTrack->GetLabel();
178     if(TMath::Abs(fLabel) > fStack->GetNtrack()) continue; 
179     
180     // read MC particle
181     AliMCParticle *mcParticle = 0x0; 
182     if(!(mcParticle = fMC->GetTrack(TMath::Abs(fLabel)))){
183       AliWarning(Form("MC particle missing for ESD fLabel %d.", fLabel));
184       continue;
185     }
186
187     AliTrackReference *ref = 0x0; 
188     Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
189     Int_t iref = 0;
190     while(iref<nRefs){
191       ref = mcParticle->GetTrackReference(iref);
192       if(ref->LocalX() > xTPC) break;
193       ref=0x0; iref++;
194     }
195
196     // read TParticle
197     TParticle *tParticle = mcParticle->Particle(); 
198     Int_t fPdg = tParticle->GetPdgCode();
199     //tParticle->IsPrimary();
200
201     //printf("[%c] ref[%2d]=", tParticle->IsPrimary() ? 'P' : 'S', iref);
202
203     TPCout=1;
204     if(ref){
205       if(ref->LocalX() > xTOF){ 
206         //printf("  TOF   [");
207         ref = mcParticle->GetTrackReference(iref-1);
208       } else {
209         //printf("%7.2f [", ref->LocalX());
210         TRDin=1;
211         if(esdTrack->GetNcls(2)) TRDout=1;
212         if(esdTrack->GetTRDpidQuality()) TRDpid=1;
213       }
214     } else { 
215       //printf("  TPC   [");
216       ref = mcParticle->GetTrackReference(iref-1);
217     }
218     Float_t pt = ref->Pt();
219     //printf("%f]\n", pt);
220     
221     TH2 *h = (TH2I*)fHistos->At(kTRDstat);
222     if(/*status & AliESDtrack::k*/TPCout) h->Fill(pt, 0);
223     if(/*status & AliESDtrack::k*/TRDin) h->Fill(pt, 1);
224     if(/*status & AliESDtrack::k*/TRDout){ 
225       ((TH1*)fHistos->At(kNCl))->Fill(esdTrack->GetNcls(2));
226       h->Fill(pt, 2);
227     }
228     if(/*status & AliESDtrack::k*/TRDpid) h->Fill(pt, 3);
229   }  
230
231   PostData(0, fHistos);
232 }
233
234
235 //____________________________________________________________________
236 void AliTRDcheckESD::Terminate(Option_t *)
237 {
238 }