fix macro control plots npart eccentricity
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTasks / AliAnalysisTaskQAflow.cxx
1 #include "TMath.h"
2 #include "TH1D.h"
3 #include "TH2D.h"
4 #include "TSeqCollection.h"
5 #include "TObjArray.h"
6 #include "TObjArray.h"
7 #include "TChain.h"
8 #include "TMCProcess.h"
9 #include "TLorentzVector.h"
10 #include "TDirectory.h"
11 #include "TROOT.h"
12 #include "TNtuple.h"
13
14 #include "AliLog.h"
15 #include "AliVParticle.h"
16 #include "AliMCParticle.h"
17 #include "AliStack.h"
18 #include "AliVEvent.h"
19 #include "AliESDEvent.h"
20 #include "AliMCEvent.h"
21 #include "AliESDtrack.h"
22 #include "AliFlowTrackCuts.h"
23 #include "AliFlowEventCuts.h"
24 #include "AliMultiplicity.h"
25 #include "AliESDtrackCuts.h"
26 #include "AliVertex.h"
27 #include "AliFlowEventSimple.h"
28 #include "AliFlowVector.h"
29
30 #include "AliAnalysisTask.h"
31 #include "AliAnalysisManager.h"
32
33 #include "AliAnalysisTaskQAflow.h"
34
35 ClassImp(AliAnalysisTaskQAflow)
36
37 //________________________________________________________________________
38 AliAnalysisTaskQAflow::AliAnalysisTaskQAflow()
39   : AliAnalysisTaskSE(),
40     fOutput(NULL),
41     fFillNtuple(kFALSE),
42     fNtuple(NULL),
43     fEventCuts(NULL),
44     fTrackCuts(NULL)
45 {
46   // Default constructor
47 }
48
49 //________________________________________________________________________
50 AliAnalysisTaskQAflow::AliAnalysisTaskQAflow(const char* name)
51   : AliAnalysisTaskSE(name),
52     fOutput(NULL),
53     fFillNtuple(kFALSE),
54     fNtuple(NULL),
55     fEventCuts(NULL),
56     fTrackCuts(NULL)
57 {
58   // Constructor
59   DefineInput(1, AliFlowEventSimple::Class());
60   DefineOutput(1, TObjArray::Class());
61 }
62
63 //________________________________________________________________________
64 void AliAnalysisTaskQAflow::UserCreateOutputObjects()
65 {
66   // Called once at the beginning
67   fOutput=new TObjArray();
68   fNtuple = new TNtuple("flowQAtree","flowQAtree","meanpt:qx:qy:mult:refmult:trig:tpcvtxx:tpcvtxy:tpcvtxz:ntrackletsA:ntrackletsC");
69   //TDirectory* before = gDirectory->mkdir("before cuts","before cuts");
70   //TDirectory* after = gDirectory->mkdir("after cuts","after cuts");
71   TObjArray* before = new TObjArray();
72   TObjArray* after = new TObjArray();
73   fOutput->Add(before);
74   fOutput->Add(after);
75
76   //define histograms
77
78   TH1* hist;  
79   hist = new TH1I("tracklet_multiplicity","tracklet multiplicity",1000,0,20000);
80   before->Add(hist); after->Add(hist->Clone()); //0
81   hist = new TH1I("track_multiplicity", "standard TPC track multiplicity",1000,0,10000);
82   before->Add(hist); after->Add(hist->Clone()); //1
83   hist = new TH1I("refmult","refmult",1000, 0,10000);
84   before->Add(hist); after->Add(hist->Clone()); //2
85   hist = new TH1I("SPD_clusters","SPD clusters",1000,0,100000);
86   before->Add(hist); after->Add(hist->Clone()); //3
87   hist = new TH1D("primary_vertexZ","primary vertex z",100,-20,20);
88   before->Add(hist); after->Add(hist->Clone()); //4
89   hist = new TH1I("ITS_clusters_on_track", "ITS clusters on track", 8, 0, 8);
90   before->Add(hist); after->Add(hist->Clone()); //5
91   hist = new TH1I("TPC_clusters_on_track", "TPC clusters on track", 159, 1, 160);
92   before->Add(hist); after->Add(hist->Clone()); //6
93   hist = new TH1D("TPC_chi2_per_cluster","TPC #chi^{2}/cluster",100,0.0,5.0);
94   before->Add(hist); after->Add(hist->Clone()); //7
95   hist = new TH1D("DCA_xy","DCA xy", 1000, -5.0, 5.0 );
96   before->Add(hist); after->Add(hist->Clone()); //8
97   hist = new TH1D("DCA_z","DCA z", 1000, -5.0, 5.0 );
98   before->Add(hist); after->Add(hist->Clone()); //9
99   hist = new TH1D("phi_tracklets","#phi tracklets", 1000, 0.0, TMath::TwoPi() );
100   before->Add(hist); after->Add(hist->Clone()); //10
101   hist = new TH1D("phi_tracks","#phi tracks", 1000, 0.0, TMath::TwoPi() );
102   before->Add(hist); after->Add(hist->Clone()); //11
103   hist = new TH1D("eta_tracklets","#eta tracklets", 1000, -2.0, 2.0 );
104   before->Add(hist); after->Add(hist->Clone()); //12
105   hist = new TH1D("eta_tracks","#eta tracks", 1000, -2.0, 2.0 );
106   before->Add(hist); after->Add(hist->Clone()); //13
107   hist = new TH1D("TPC_vertex_z", "TPC vertex z", 100,-20.0,20.0);
108   before->Add(hist); after->Add(hist->Clone()); //14
109   hist = new TH1D("ptyield", "p_{t} spectrum", 10000,0.0,10.0);
110   before->Add(hist); after->Add(hist->Clone()); //15
111   
112   //post data here as it doesn't change anyway (the pointer to list anyway)
113
114   //restore dir add status
115   PostData(1, fOutput);
116 }
117
118 //________________________________________________________________________
119 void AliAnalysisTaskQAflow::UserExec(Option_t *)
120 {
121
122   //get teh input data
123   AliESDEvent* event = dynamic_cast<AliESDEvent*>(InputEvent());
124   if (!event)
125   {
126     AliFatal("no ESD event");
127     return;
128   }
129
130   //TObjArray* before = ((TDirectory*)fOutput->At(0))->GetList();
131   //TObjArray* after = ((TDirectory*)fOutput->At(1))->GetList();
132   TObjArray* before = (TObjArray*)fOutput->At(0);
133   TObjArray* after = (TObjArray*)fOutput->At(1);
134   TH1* htrackletmultB = dynamic_cast<TH1*>(before->At(0));
135   TH1* htrackletmultA = dynamic_cast<TH1*>(after->At(0));
136   TH1* htrackmultB = dynamic_cast<TH1*>(before->At(1));
137   TH1* htrackmultA = dynamic_cast<TH1*>(after->At(1));
138   TH1* hrefmultB = dynamic_cast<TH1*>(before->At(2));
139   TH1* hrefmultA = dynamic_cast<TH1*>(after->At(2));
140   TH1* hspdclustersB = dynamic_cast<TH1*>(before->At(3));
141   TH1* hspdclustersA = dynamic_cast<TH1*>(after->At(3));
142   TH1* hprimvtxzB = dynamic_cast<TH1*>(before->At(4));
143   TH1* hprimvtxzA = dynamic_cast<TH1*>(after->At(4));
144
145   TH1* hITSclsB = dynamic_cast<TH1*>(before->At(5));
146   TH1* hITSclsA = dynamic_cast<TH1*>(after->At(5));
147   TH1* hTPCclsB = dynamic_cast<TH1*>(before->At(6));
148   TH1* hTPCclsA = dynamic_cast<TH1*>(after->At(6));
149   TH1* hTPCchi2B = dynamic_cast<TH1*>(before->At(7));
150   TH1* hTPCchi2A = dynamic_cast<TH1*>(after->At(7));
151   TH1* hdcaxyB = dynamic_cast<TH1*>(before->At(8));
152   TH1* hdcaxyA = dynamic_cast<TH1*>(after->At(8));
153   TH1* hdcazB = dynamic_cast<TH1*>(before->At(9));
154   TH1* hdcazA = dynamic_cast<TH1*>(after->At(9));
155   TH1* hphitrackletsB = dynamic_cast<TH1*>(before->At(10));
156   TH1* hphitrackletsA = dynamic_cast<TH1*>(after->At(10));
157   TH1* hphitracksB = dynamic_cast<TH1*>(before->At(11));
158   TH1* hphitracksA = dynamic_cast<TH1*>(after->At(11));
159   TH1* hetatrackletsB = dynamic_cast<TH1*>(before->At(12));
160   TH1* hetatrackletsA = dynamic_cast<TH1*>(after->At(12));
161   TH1* hetatracksB = dynamic_cast<TH1*>(before->At(13));
162   TH1* hetatracksA = dynamic_cast<TH1*>(after->At(13));
163   TH1* hprimvtxzTPCB = dynamic_cast<TH1*>(before->At(14));
164   TH1* hprimvtxzTPCA = dynamic_cast<TH1*>(after->At(14));
165   TH1* hptyieldB = dynamic_cast<TH1*>(before->At(15));
166   TH1* hptyieldA = dynamic_cast<TH1*>(after->At(15));
167
168   AliMultiplicity* tracklets = const_cast<AliMultiplicity*>(event->GetMultiplicity());
169   Int_t ntracklets=0;
170   Int_t nspdclusters=0;
171   Int_t nspd1clusters=0;
172   Int_t ntrackletsA=0;
173   Int_t ntrackletsC=0;
174   if (tracklets)
175   {
176     ntracklets = tracklets->GetNumberOfTracklets();
177     nspdclusters = tracklets->GetNumberOfITSClusters(0,1);
178     nspd1clusters = tracklets->GetNumberOfITSClusters(1);
179     for (Int_t i=0; i<tracklets->GetNumberOfTracklets(); i++)
180     {
181       Bool_t pass=fTrackCuts->IsSelected(tracklets,i);
182       Float_t phi=tracklets->GetPhi(i);
183       Float_t eta=tracklets->GetEta(i);
184       hphitrackletsB->Fill(phi); if (pass) hphitrackletsA->Fill(phi);
185       hetatrackletsB->Fill(eta); if (pass) hetatrackletsA->Fill(eta); 
186       if (eta>0) ntrackletsC++;
187       else ntrackletsA++;
188     }
189     
190   }
191   //Int_t trackmult = AliESDtrackCuts::GetReferenceMultiplicity(event,kTRUE);
192   AliFlowEventCuts newcuts;
193   newcuts.SetRefMultCuts(fTrackCuts);
194   Int_t trackmult = newcuts.RefMult(event);
195   Int_t refmult = fEventCuts->RefMult(event);
196   Bool_t passevent = fEventCuts->IsSelected(event);
197   htrackletmultB->Fill(ntracklets); if (passevent) htrackletmultA->Fill(ntracklets); 
198   htrackmultB->Fill(trackmult); if (passevent) htrackmultA->Fill( trackmult); 
199   hrefmultB->Fill(refmult); if (passevent) hrefmultA->Fill( refmult); 
200   hspdclustersB->Fill(nspdclusters); if (passevent) hspdclustersA->Fill( nspdclusters);
201   const AliVertex* vertex = event->GetPrimaryVertex();
202   Float_t vtxz=0.0;
203   Float_t vtxx=0.0;
204   Float_t vtxy=0.0;
205   if (vertex)
206   {
207     vtxz = vertex->GetZ();
208     vtxx = vertex->GetX();
209     vtxy = vertex->GetY();
210     hprimvtxzB->Fill(vtxz); if (passevent) hprimvtxzA->Fill(vtxz);
211   }
212   const AliVertex* vertextpc = event->GetPrimaryVertexTPC();
213   Float_t vtxTPCx=0.0;
214   Float_t vtxTPCy=0.0;
215   Float_t vtxTPCz=0.0;
216   if (vertextpc)
217   {
218     vtxTPCx = vertextpc->GetX();
219     vtxTPCy = vertextpc->GetY();
220     vtxTPCz = vertextpc->GetZ();
221     hprimvtxzTPCB->Fill(vtxTPCz); if (passevent) hprimvtxzTPCA->Fill(vtxTPCz);
222   }
223   fTrackCuts->SetEvent(event);
224   Float_t meanpt=0;
225   Int_t ntracks=fTrackCuts->GetNumberOfInputObjects();
226   for (Int_t i=0; i<ntracks; i++)
227   {
228     TObject* obj = fTrackCuts->GetInputObject(i);
229     Bool_t pass = fTrackCuts->IsSelected(obj,i);
230     Float_t dcaxy=0.0;
231     Float_t dcaz=0.0;
232     Float_t tpcchi2=0.0;
233     Float_t tpcchi2percls=0.0;
234     Int_t ntpccls=0;
235     Int_t nitscls=0;
236     Float_t eta=0.0;
237     Float_t phi=0.0;
238     Float_t pt=0.0;
239     AliESDtrack* track = dynamic_cast<AliESDtrack*>(fTrackCuts->GetTrack());
240     if (track)
241     {
242       track->GetImpactParameters(dcaxy,dcaz);
243       tpcchi2=track->GetTPCchi2();
244       ntpccls=track->GetTPCNcls();
245       eta=track->Eta();
246       phi=track->Phi();
247       pt=track->Pt();
248       meanpt+=pt;
249       tpcchi2percls= (ntpccls==0)?0.0:tpcchi2/ntpccls;
250       nitscls = track->GetNcls(0);
251       hITSclsB->Fill(nitscls); if (pass) hITSclsA->Fill( nitscls);
252       hTPCclsB->Fill(ntpccls); if (pass) hTPCclsA->Fill( ntpccls);
253       hTPCchi2B->Fill(tpcchi2percls); if (pass) hTPCchi2A->Fill( tpcchi2percls);
254       hdcaxyB->Fill(dcaxy); if (pass) hdcaxyA->Fill( dcaxy);
255       hdcazB->Fill(dcaz); if (pass) hdcazA->Fill(dcaz);
256       hetatracksB->Fill(eta); if (pass) hetatracksA->Fill(eta);
257       hphitracksB->Fill(phi); if (pass) hphitracksA->Fill(phi);
258       hptyieldB->Fill(pt); if (pass) hptyieldA->Fill(pt);
259     }
260   }
261   meanpt = meanpt/ntracks;
262
263   ///////////////////////////////////////////////////////////////////////
264   //the ntuple part/////////////
265   AliFlowEventSimple* flowevent = dynamic_cast<AliFlowEventSimple*>(GetInputData(1));
266   if (flowevent && fFillNtuple)
267   {
268     AliFlowVector qvec = flowevent->GetQ(2);
269     Double_t qx = qvec.X();
270     Double_t qy = qvec.Y();
271     TString triggers = event->GetFiredTriggerClasses();
272     Int_t trig=0;
273     if (triggers.Contains("CMBAC-B-NOPF-ALL")) trig+=1;
274     if (triggers.Contains("CMBS2C-B-NOPF-ALL")) trig+=10;
275     if (triggers.Contains("CMBS2A-B-NOPF-ALL")) trig+=100;
276     fNtuple->Fill(meanpt,qx,qy,trackmult,refmult,trig,vtxTPCx,vtxTPCy,vtxTPCz,ntrackletsA,ntrackletsC);
277   }//if flowevent
278 }
279
280 //________________________________________________________________________
281 void AliAnalysisTaskQAflow::Terminate(Option_t *)
282 {
283   
284 }
285
286 //________________________________________________________________________
287 AliAnalysisTaskQAflow::~AliAnalysisTaskQAflow()
288 {
289   //dtor
290   delete fTrackCuts;
291   delete fEventCuts;
292 }