Adding analysis for Flat ESD case
[u/mrichter/AliRoot.git] / TPC / Calib / AliAnalysisTaskPt.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TH1F.h"
4 #include "TCanvas.h"
5 #include "TObjArray.h"
6
7 #include "AliAnalysisTask.h"
8 #include "AliAnalysisManager.h"
9
10 #include "AliESDEvent.h"
11 #include "AliVVevent.h"
12 #include "AliVVtrack.h"
13 #include "AliESDtrackCuts.h"
14 #include "AliVEventHandler.h"
15 #include "AliTPCseed.h"
16 #include "AliTPCclusterMI.h"
17 #include "AliVVfriendEvent.h"
18 #include "AliVVfriendTrack.h"
19
20 #include "AliAnalysisTaskPt.h"
21
22 // example of an analysis task creating a p_t spectrum
23 // Authors: Panos Cristakoglou, Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing
24
25 ClassImp(AliAnalysisTaskPt)
26
27 //________________________________________________________________________
28 AliAnalysisTaskPt::AliAnalysisTaskPt(const char *name) 
29 : AliAnalysisTask(name, ""), fESD(0), fESDfriend(0), fHistPt(0), fCuts(0), fEv(0), fHistQ(0), fListOut(0), fUseFriends(kFALSE), fHistNTPCCl(0), fHistNESDtracks(0),   fHistNESDfriendtracks(0)
30
31 {
32   // Constructor
33
34   // Define input and output slots here
35   // Input slot #0 works with a TChain
36   DefineInput(0, TChain::Class());
37   // Output slot #0 writes into a TH1 container
38   DefineOutput(0, TList::Class());
39 }
40
41 //________________________________________________________________________
42 void AliAnalysisTaskPt::ConnectInputData(Option_t *) 
43 {
44   // Connect ESD or AOD here
45   // Called once
46
47   printf("----> AliAnalysisTaskPt::ConnectInputData\n");
48   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
49   if (!tree) {
50     Printf("ERROR: Could not read chain from input slot 0");
51   } else {
52     // Disable all branches and enable only the needed ones
53     // The next two lines are different when data produced as AliESDEvent is read
54     /*
55     tree->SetBranchStatus("*", kFALSE);
56     tree->SetBranchStatus("fTracks.*", kTRUE);
57     */
58
59     AliVEventHandler *esdH = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
60
61     if (!esdH) {
62       Printf("ERROR: Could not get ESDInputHandler");
63     } else {
64       Printf("----> AliAnalysisTaskPt::ConnectInputData Getting the event from handler %p", esdH);
65       //fESD = dynamic_cast<AliESDEvent*>(esdH->GetEvent());
66       fESD = esdH->GetVVEvent();
67       if (fUseFriends){ 
68         fESDfriend = esdH->GetVVFriendEvent();
69       }
70     }
71     if (!fESD) {
72       Printf("ERROR, no ESD event");
73     }
74     if (fUseFriends && !fESDfriend) {
75       Printf("ERROR, no ESD friend");
76     }
77   }
78
79   Printf("fESD = %p, fESDfriend = %p", fESD, fESDfriend);
80   printf("<---- AliAnalysisTaskPt::ConnectInputData\n");
81 }
82
83 //________________________________________________________________________
84 void AliAnalysisTaskPt::CreateOutputObjects()
85 {
86   // Create histograms
87   // Called once
88
89   fListOut = new TList();
90   fListOut->SetOwner();
91   fListOut->SetName("listHistos");
92
93   fHistPt = new TH1F("fHistPt", "P_{T} distribution", 15, 0.1, 3.1);
94   fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
95   fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
96   fHistPt->SetMarkerStyle(kFullCircle);
97
98   fHistQ = new TH1F("fHistQ", "TPC clusters: Q distribution", 1000, 0, 10000);
99   fHistQ->GetXaxis()->SetTitle("Q");
100   fHistQ->GetYaxis()->SetTitle("dN/dQ");
101   fHistQ->SetMarkerStyle(kFullCircle);
102
103   fHistNTPCCl = new TH1F("fHistNTPCCl", "Number of TPC clusters", 160, -0.5, 159.5);
104   fHistNTPCCl->GetXaxis()->SetTitle("n. TPC Cl.");
105   fHistNTPCCl->GetYaxis()->SetTitle("dN/d(n. TPC Cl)");
106   fHistNTPCCl->SetMarkerStyle(kFullCircle);
107
108   fHistNESDtracks = new TH1F("fHistNESDtracks", "Number of ESD friend tracks", 1000, -0.5, 999.5);
109   fHistNESDtracks->GetXaxis()->SetTitle("n. ESD friend tracks");
110   fHistNESDtracks->GetYaxis()->SetTitle("dN/d(n. ESD friend tracks)");
111   fHistNESDtracks->SetMarkerStyle(kFullCircle);
112
113   fHistNESDfriendtracks = new TH1F("fHistNESDfriendtracks", "Number of ESD tracks", 1000, -0.5, 999.5);
114   fHistNESDfriendtracks->GetXaxis()->SetTitle("n. ESD tracks");
115   fHistNESDfriendtracks->GetYaxis()->SetTitle("dN/d(n. ESD tracks)");
116   fHistNESDfriendtracks->SetMarkerStyle(kFullCircle);
117
118   fListOut->Add(fHistPt);
119   fListOut->Add(fHistQ);
120   fListOut->Add(fHistNTPCCl);
121   fListOut->Add(fHistNESDtracks);
122   fListOut->Add(fHistNESDfriendtracks);
123
124   PostData(0, fListOut);
125
126   fCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(1);
127 }
128
129 //________________________________________________________________________
130 void AliAnalysisTaskPt::Exec(Option_t *) 
131 {
132   // Main loop
133   // Called for each event
134
135   if (!fESD) {
136     Printf("ERROR: fESD not available");
137     return;
138   }
139   if (!fESDfriend) {
140     Printf("ERROR: fESDfriend not available");
141       if (fUseFriends){
142         return;
143       }
144   }
145
146   Int_t nESDtracks = fESD->GetNumberOfTracks();
147   Int_t nESDfriendtracks = 0;
148   if (fUseFriends) nESDfriendtracks = fESDfriend->GetNumberOfTracks();
149   Printf("There are %d tracks in this event", nESDtracks);
150   Printf("... and there are %d friends in this event", nESDfriendtracks);
151
152   fHistNESDtracks->Fill(nESDtracks);
153   fHistNESDfriendtracks->Fill(nESDfriendtracks);
154
155   // Track loop to fill a pT spectrum
156   for (Int_t iTracks = 0; iTracks < nESDtracks; iTracks++) {
157     Printf("Checking track %d: Note that with Flat, the GetTrack is not yet implemented!!!", iTracks);
158     AliVVtrack* track = fESD->GetVVTrack(iTracks);
159     if (!track) {
160       Printf("ERROR: Could not receive track %d", iTracks);
161       continue;
162     }
163     Printf("track %d has pt = %f", iTracks, track->GetPt());
164     fHistPt->Fill(track->GetPt());
165     fHistNTPCCl->Fill(track->GetTPCNcls());
166   } //track loop 
167
168
169   if (fUseFriends){
170     // Friend Track loop
171     for (Int_t iFriend = 0; iFriend < nESDfriendtracks; iFriend++) {
172       //Printf("Getting friend %d", iFriend);
173       const AliVVfriendTrack* friendTrack = fESDfriend->GetTrack(iFriend);
174       if (!friendTrack) {
175         Printf("ERROR: Could not receive track %d", iFriend);
176         continue;
177       }
178       
179       AliTPCseed seed;
180       Int_t err = friendTrack->GetTPCseed( seed );
181       if( err==0 ){
182         Printf("Found TPC seed" );
183         for (Int_t irow = 0; irow < 160; irow++){
184           AliTPCclusterMI* cluMI = seed.GetClusterPointer(irow);
185           if (cluMI){
186             Printf("Found cluster at row %d", irow);
187             Float_t q = cluMI->GetQ();
188             Printf("Q = %f", q);
189             fHistQ->Fill(q);
190           }
191           else {
192             Printf("Row %d does not contain clusters", irow);
193           }              
194         }
195       }     
196       else {
197         //Printf("Schade... seed is %p", seed);
198       }
199     }
200   }
201
202   // Post output data.
203   PostData(0, fListOut);
204   fEv++;
205 }      
206
207 //________________________________________________________________________
208 void AliAnalysisTaskPt::Terminate(Option_t *) 
209 {
210   // Draw result to the screen
211   // Called once at the end of the query
212
213   Printf("Terminate called: fESD = %p", fESD);
214
215   fListOut = dynamic_cast<TList*> (GetOutputData(0)); 
216
217   if (fListOut) {
218     fHistPt = dynamic_cast<TH1F*>(fListOut->FindObject("fHistPt")); 
219     if (!fHistPt) {
220       Printf("ERROR: fHistPt not available");
221       return;
222     }
223    
224     TCanvas *c1 = new TCanvas("AliAnalysisTaskPt","Pt",10,10,510,510);
225     c1->cd(1)->SetLogy();
226     fHistPt->DrawCopy("E");
227   }
228   else {
229     Printf("In Terminate: no TList found");
230   }
231
232 }