In AliMUONClusterInfo:
[u/mrichter/AliRoot.git] / PWGUD / FirstPhysics / AliAnalysisTaskFPexample.cxx
1
2 #include "AliAnalysisTaskFPexample.h"
3
4 #include "TH1D.h"
5 #include "TH2D.h"
6 #include "TCanvas.h"
7 #include "TList.h"
8 #include "TString.h"
9
10 #include "AliAnalysisManager.h"
11 #include "AliESDtrackCuts.h"
12 #include "AliESDEvent.h"
13 #include "AliESDHeader.h"
14 #include "AliMCEvent.h"
15 #include "AliMultiplicity.h"
16 #include "AliTriggerAnalysis.h"
17
18 ClassImp(AliAnalysisTaskFPexample)
19
20 //________________________________________________________________________
21 AliAnalysisTaskFPexample::AliAnalysisTaskFPexample(const char *name) :
22   AliAnalysisTaskFirstPhysics(name),
23   fhTrackPt(0),
24   fh2TrackPhiEta(0),
25   fhMulITSTPC(0),
26   fhMulITSSA(0),
27   fhMulSPD(0),
28   fh2TrackletsPhiEta(0),
29   fh2TracksPhiTPCchi2(0)
30 {
31   DefineOutput(1, TList::Class());
32 }
33
34 //________________________________________________________________________
35 AliAnalysisTaskFPexample::~AliAnalysisTaskFPexample()
36 {
37   // Destructor. Clean-up the output list, but not the histograms that are put inside
38   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
39   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
40     delete fOutput;
41   }
42   for (Int_t i = 0; i < knTrackCuts; i ++) {
43     delete fTrackCuts[i];
44     fTrackCuts[i] = 0;
45   }
46 }
47
48 //________________________________________________________________________
49 void AliAnalysisTaskFPexample::UserCreateOutputObjects()
50 {
51   PrepareOutputList();
52   // Define cuts
53   PrepareDefaultTrackCuts();
54
55   // Create histograms
56   const Int_t ptbins = 15;
57   const Double_t ptlow = 0.1, ptup = 3.1;
58   const Int_t etabins = 40;
59   const Double_t etalow = -2.0, etaup = 2.0;
60   const Int_t mulbins = 200;
61   const Int_t phibins = 120;
62
63   fhTrackPt = UserHisto1d("fhTrackPt", "p_{T} distribution for ESD",
64                           "P_{T} (GeV/c)", ptbins, ptlow, ptup);
65   fh2TrackPhiEta = UserHisto2d("fh2TrackPhiEta", "ESD tracks",
66                                "#Phi", phibins, 0, 2 * TMath::Pi(),
67                                "#eta", etabins, etalow, etaup);
68   fhMulITSTPC = UserHisto1d("fhMulITSTPC", "N_{CH} distribution for ITS+TPC tracks",
69                             "N_{CH}", mulbins, 0, mulbins);
70   fhMulITSSA = UserHisto1d("fhMulITSSA", "N_{CH} distribution for ITS SA tracks",
71                            "N_{CH}", mulbins, 0, mulbins);
72   fhMulSPD = UserHisto1d("fhMulSPD", "N_{CH} distribution for SPD tracklets",
73                          "N_{CH}", mulbins, 0, mulbins);
74   fh2TrackletsPhiEta = UserHisto2d("fh2TrackletsPhiEta", "Tracklets",
75                                    "#Phi", phibins, 0, 2 * TMath::Pi(),
76                                    "#eta", etabins, etalow, etaup);
77   fh2TracksPhiTPCchi2 = UserHisto2d("fh2TracksPhiTPCchi2", "ESD tracks #Phi vs. TPC #chi^{2}",
78                                     "#Phi", phibins, 0, 2 * TMath::Pi(),
79                                     "#chi^{2}", 500, 0, 500);
80
81   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
82 }
83
84 //________________________________________________________________________
85 void AliAnalysisTaskFPexample::UserExec(Option_t *)
86 {
87   if (!GetESDEvent()) {
88     AliError("Unable to read the ESD");
89     return;
90   }
91   if (!CheckVertex()) {
92     return;
93   }
94
95   const Int_t nESDTracks = fESD->GetNumberOfTracks();
96
97   Int_t iHighestID = 0;
98   for (Int_t iTrack = 0; iTrack < nESDTracks; iTrack++) {
99     if (fESD->GetTrack(iTrack)->GetLabel() > iHighestID) {
100       iHighestID = fESD->GetTrack(iTrack)->GetLabel();
101     }
102   }
103   const Int_t nMaxID = iHighestID + 1;
104   bool aGlobalBits[nMaxID], aPureITSBits[nMaxID];
105   for (Int_t iFlag = 0; iFlag < nMaxID; iFlag++) {
106     aGlobalBits[iFlag] = false;
107     aPureITSBits[iFlag] = false;
108   }
109
110   // flags for secondary and rejected tracks
111   const int kRejectBit = BIT(15); // set this bit in ESD tracks if it is rejected by a cut
112   const int kSecondaryBit = BIT(16); // set this bit in ESD tracks if it is secondary according to a cut
113
114   Int_t nTracksITSTPC = 0; // multiplicity counters
115   Int_t nTracksITSSA = 0;
116   Int_t nTrackletsSPD = 0;
117
118   for(Int_t iTrack = 0; iTrack < nESDTracks; iTrack++){
119     AliESDtrack* esdTrack = fESD->GetTrack(iTrack);
120     // if track is a secondary from a V0, flag as a secondary
121     if (esdTrack->IsOn(AliESDtrack::kMultInV0)) {
122       esdTrack->SetBit(kSecondaryBit);
123       continue;
124     }
125     AliESDtrack* tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(fESD, esdTrack->GetID());
126     // check tracks with ITS part
127     if (esdTrack->IsOn(AliESDtrack::kITSin)) {
128       if (!esdTrack->IsOn(AliESDtrack::kITSpureSA)) { // track has ITS part but is not an ITS_SA
129         // TPC+ITS
130         if (esdTrack->IsOn(AliESDtrack::kTPCin)) {  // Global track, has ITS and TPC contributions
131           if (fTrackCuts[kTrackCutQGlo]->AcceptTrack(esdTrack)) { // good ITS+TPC track
132             if (fTrackCuts[kTrackCutDCAwSPD]->AcceptTrack(esdTrack)
133                 || fTrackCuts[kTrackCutDCAwoSPD]->AcceptTrack(esdTrack)) {
134               nTracksITSTPC++;
135               fhTrackPt->Fill(esdTrack->Pt());
136               fh2TrackPhiEta->Fill(esdTrack->Phi(), esdTrack->Eta());
137               if (tpcTrack) {
138                 fh2TracksPhiTPCchi2->Fill(esdTrack->Phi(), tpcTrack->GetTPCchi2());
139               }
140             } else {
141               esdTrack->SetBit(kSecondaryBit); // large DCA -> secondary, don't count either track not associated tracklet
142             }
143           } else {
144             esdTrack->SetBit(kRejectBit); // bad quality, don't count the track, but may count tracklet if associated
145           }
146         } else if (fTrackCuts[kTrackCutQITS]->AcceptTrack(esdTrack)) { // good ITS complementary track
147           if (fTrackCuts[kTrackCutDCAwSPD]->AcceptTrack(esdTrack)
148               || fTrackCuts[kTrackCutDCAwoSPD]->AcceptTrack(esdTrack)) {
149             nTracksITSTPC++;
150             fhTrackPt->Fill(esdTrack->Pt());
151             fh2TrackPhiEta->Fill(esdTrack->Phi(), esdTrack->Eta());
152             if (tpcTrack) {
153               fh2TracksPhiTPCchi2->Fill(esdTrack->Phi(), tpcTrack->GetTPCchi2());
154             }
155           } else {
156             esdTrack->SetBit(kSecondaryBit); // large DCA -> secondary, don't count either track not associated tracklet
157           }
158         } else {
159           esdTrack->SetBit(kRejectBit); // bad quality, don't count the track, but may count tracklet if associated
160         }
161       } else { // pure ITS SA tracks
162         if (fTrackCuts[kTrackCutQITS]->AcceptTrack(esdTrack)) { // good ITSSA track
163           if (fTrackCuts[kTrackCutDCAwSPD]->AcceptTrack(esdTrack)
164               || fTrackCuts[kTrackCutDCAwoSPD]->AcceptTrack(esdTrack)) {
165             nTracksITSSA++;
166           } else {
167             esdTrack->SetBit(kRejectBit);
168           }
169         } else {
170           esdTrack->SetBit(kRejectBit);
171         }
172       }
173     }
174     if (tpcTrack) {
175       delete tpcTrack;
176     }
177   }
178
179   // get multiplicity from ITS tracklets to complement TPC+ITS, and ITS pure SA
180   const AliMultiplicity* multiplicitySPD = fESD->GetMultiplicity();    // spd multiplicity object
181   Int_t id1, id2, id3, id4;
182   AliESDtrack *tr1 = 0, *tr3 = 0;
183   for (Int_t iTracklet = 0; iTracklet < multiplicitySPD->GetNumberOfTracklets(); iTracklet++) {
184     if (TMath::Abs(multiplicitySPD->GetEta(iTracklet)) > GetCutEta()) {
185       continue; // eta selection for tracklets
186     }
187     nTrackletsSPD++;
188     fh2TrackletsPhiEta->Fill(multiplicitySPD->GetPhi(iTracklet), multiplicitySPD->GetEta(iTracklet));
189     // if counting tracks+tracklets, check if clusters were already used in tracks
190     // and get the id of the tracks in which they were used
191
192     // references for eventual Global/ITS_SA tracks
193     multiplicitySPD->GetTrackletTrackIDs(iTracklet, 0, id1, id2);
194     tr1 = id1 >= 0 ? fESD->GetTrack(id1) : 0;
195
196     // references for eventual ITS_SA_pure tracks
197     multiplicitySPD->GetTrackletTrackIDs(iTracklet, 1, id3, id4);
198     tr3 = id3 >= 0 ? fESD->GetTrack(id3) : 0;
199
200     // are both clusters from the same tracks? If not, skip the
201     // tracklet (shouldn't change things much)
202     if (id1 != id2 || id3 != id4) {
203       continue;
204     }
205
206     // has associated global track been associated to a previous tracklet?
207     bool bUsedInGlobal = (id1 != -1) ? aGlobalBits[id1] : false;
208     // has associated pure ITS track been associated to a previous tracklet?
209     bool bUsedInPureITS = (id3 != -1) ? aPureITSBits[id3] : false;
210
211     // counting tracklet as global+complementary track
212     if ((tr1 && !tr1->TestBit(kSecondaryBit)) // reject as secondary
213         && (tr1 && tr1->TestBit(kRejectBit))) { // already accounted
214       if (!bUsedInGlobal) {
215         nTracksITSTPC++;
216         if (id1 > 0) {
217           // mark global track linked to this tracklet as associated
218           aGlobalBits[id1] = true;
219         }
220       }
221     } else if (id1 < 0) {
222       nTracksITSTPC++;
223     }
224
225     // counting tracklet as ITS SA pure track
226     if ((tr3 && tr3->TestBit(kSecondaryBit))
227         && (tr3 && !tr3->TestBit(kRejectBit))) {
228       if (!bUsedInPureITS) {
229         nTracksITSSA++;
230         if (id3 > 0) {
231           // mark global track linked to this tracklet as associated
232           aPureITSBits[id3] = true;
233         }
234       }
235     } else if (id3 < 0) {
236       nTracksITSSA++;
237     }
238
239   }
240   fhMulITSTPC->Fill(nTracksITSTPC);
241   fhMulITSSA->Fill(nTracksITSSA);
242   fhMulSPD->Fill(nTrackletsSPD);
243
244   PostData(1, fOutput);
245 }
246
247 //________________________________________________________________________
248 void AliAnalysisTaskFPexample::Terminate(Option_t *)
249 {
250   // Draw result to screen, or perform fitting, normalizations
251   // don't get too fancy, keep your histos raw and massage them with macros
252
253   fOutput = dynamic_cast<TList*> (GetOutputData(1));
254   if(!fOutput) {
255     Printf("ERROR: could not retrieve TList fOutput");
256     return;
257   }
258
259   if (!(GetHisto1FromOutput("fhTrackPt", fhTrackPt) &&
260         GetHisto2FromOutput("fh2TrackPhiEta", fh2TrackPhiEta) &&
261         GetHisto1FromOutput("fhMulITSTPC", fhMulITSTPC) &&
262         GetHisto1FromOutput("fhMulITSSA", fhMulITSSA) &&
263         GetHisto1FromOutput("fhMulSPD", fhMulSPD) &&
264         GetHisto2FromOutput("fh2TrackletsPhiEta", fh2TrackletsPhiEta) &&
265         GetHisto2FromOutput("fh2TracksPhiTPCchi2", fh2TracksPhiTPCchi2))) {
266     AliError("Couldn't load every histogram from output.");
267     return;
268   }
269
270   TCanvas *c = new TCanvas("AliAnalysisTaskFPexample", "Data Quality Quick Overview");
271   c->Divide(2, 2);
272   c->cd(1)->SetLogy();
273   fhTrackPt->DrawCopy("E");
274   c->cd(2);
275   fh2TrackPhiEta->DrawCopy("");
276   c->cd(3)->SetLogy();
277   fhMulITSTPC->DrawCopy("E");
278   c->cd(4)->SetLogy();
279   fhMulITSSA->DrawCopy("");
280 }