Removing useless flag.
[u/mrichter/AliRoot.git] / T0 / AliT0MultiplicityTask.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TH1F.h"
4 #include "TH2F.h"
5 #include "TList.h"
6 #include "TClonesArray.h"
7 #include "TCanvas.h"
8 #include "TStyle.h"
9 #include "TLine.h"
10 #include "TText.h"
11 #include "TFile.h"
12 #include "TBenchmark.h"
13
14 #include "AliAnalysisTask.h"
15 #include "AliAnalysisManager.h"
16
17 #include "AliLog.h"
18 #include "AliVEvent.h"
19 #include "AliESDEvent.h"
20 #include "AliMCEvent.h"
21 #include "AliESDtrack.h"
22 #include "AliStack.h"
23 #include "AliTrackReference.h"
24
25 #include "AliT0MultiplicityTask.h"
26
27 ClassImp(AliT0MultiplicityTask)
28
29
30 AliT0MultiplicityTask::AliT0MultiplicityTask()
31   : AliAnalysisTaskSE("AliT0MultiplicityTask"),
32     fListOfHistos(0),
33     fOrA(0),
34     fOrC(0),
35     fMean(0),
36     fVertex(0),
37     fTime(0),
38     fAmp(0),
39     fTotalMult(0),
40     fMultRecTotal(0),
41     fMultRecRealA(0),
42     fMultRecRealC(0),
43     fPrim(0),
44     fRef(0),
45     fRec(0)
46 {
47   // Default constructor
48   AliInfo("Default constructor AliT0MultiplicityTask");
49   // Define input and output slots here
50   // Input slot #0 works with a TChain
51   DefineInput(0, TChain::Class());
52   // Output slot #1 TList
53   DefineOutput(1, TList::Class());
54 }
55
56
57 AliT0MultiplicityTask::AliT0MultiplicityTask(const char* name)
58   : AliAnalysisTaskSE(name),
59     fListOfHistos(0),
60     fOrA(0),
61     fOrC(0),
62     fMean(0),
63     fVertex(0),
64     fTime(0),
65     fAmp(0),   
66     fTotalMult(0),
67     fMultRecTotal(0),
68     fMultRecRealA(0),
69     fMultRecRealC(0),
70     fPrim(0),
71     fRef(0),
72     fRec(0)
73 {
74   // Constructor
75   // Define input and output slots here
76   // Input slot #0 works with a TChain
77   DefineInput(0, TChain::Class());
78   // Output slot #1 TList
79   DefineOutput(1, TList::Class());
80 }
81
82
83 void AliT0MultiplicityTask::UserCreateOutputObjects() 
84 {
85   // Create histograms
86   // Called once
87   // Create output container
88   fListOfHistos = new TList();
89   
90  
91   fOrA = new TH1F("hOrA"," T0 A", 100,2450,2700);
92   fOrC = new TH1F("hOrC"," T0 C", 100,2450,2700);
93   fMean= new TH1F("hMean"," T0 ",100, 12000, 13000);
94   fVertex = new TH1F("hVertex","Z position of vertex",   100,-30, 30);
95    
96   fTime = new TH1F("Time", " Amp vs Time",100, 12000, 13000);
97   fAmp = new TH1F("fAmp", " Amplitude", 100, 0, 200);
98
99   //  fHighMult  = new TH1F("fHighMult", " events with high amp ",100, 0, 10);
100   
101   fTotalMult = new TH1F("fTotalMult","total multiplicity",500,0,5000);
102   
103   fMultRecRealA = new TH2F("fMultRecRealA"," ",100,0.,200,100,0,200); 
104   fMultRecRealC = new TH2F("fMultRecRealC"," ",100,0,200,100,0,200); 
105   fMultRecTotal = new TH2F("fMultRecTotal"," ",100,0,200,100,0,200); 
106   
107   fPrim = new TH1F("fPrim", " primary",100, 0, 200);
108   fRef  = new TH1F("fRef", " from TR ",100, 0, 200);
109   fRec  = new TH1F("fRec", " in ESD ",100, 0, 200);
110    
111   fListOfHistos->Add(fOrA);
112   fListOfHistos->Add(fOrC);
113   fListOfHistos->Add(fMean);
114   fListOfHistos->Add(fVertex);
115   fListOfHistos->Add(fAmp);
116   fListOfHistos->Add(fTime);
117   fListOfHistos->Add(fTotalMult);
118   fListOfHistos->Add(fMultRecRealA);
119   fListOfHistos->Add(fMultRecRealC);
120   fListOfHistos->Add(fMultRecTotal);
121   fListOfHistos->Add(fPrim);
122   fListOfHistos->Add(fRef);
123   fListOfHistos->Add(fRec);
124 }  
125
126
127 void AliT0MultiplicityTask::UserExec(Option_t *) 
128 {
129   // Main loop
130   // Called for each event
131
132     
133   // MC information
134   AliMCEvent* mcEvent = MCEvent();
135   if (!mcEvent) {
136     Printf("ERROR: Could not retrieve MC event");
137     return;
138   }
139   
140   
141   Int_t primaryA=0;
142   Int_t primaryC=0;
143   Int_t numPrim = 0;
144   Int_t refT0=0;
145
146   AliStack* stack = mcEvent->Stack();
147   // printf("AliT0MultiplicityAnalysis: Number of tracks on stack %5d\n", stack->GetNtrack());
148   Int_t nTracks  = mcEvent->GetNumberOfTracks();
149   for (Int_t ipart=0; ipart<nTracks; ipart++)
150     {
151       //    TParticle* particle = stack->Particle(ipart);
152       AliMCParticle* track = mcEvent->GetTrack(ipart);
153       if (!track) {
154         Printf("ERROR: Could not receive track %d (mc loop)", ipart);
155         continue;
156       }
157       Int_t label = track->GetLabel();
158       if(stack->IsPhysicalPrimary(label) == kFALSE)
159         continue;
160       
161       numPrim++;
162       Double_t eta=track->Eta();
163       
164       if(eta<-2.97 && eta>-3.28)  primaryC++;
165       if (eta >4.61 && eta<4.92)  primaryA++;
166     
167        // Loop over Track References
168        AliTrackReference* trackRef = 0;
169
170       for (Int_t iTrackRef = 0; iTrackRef  < track->GetNumberOfTrackReferences(); iTrackRef++) {
171         trackRef = track->GetTrackReference(iTrackRef);
172         if(trackRef) {
173           Int_t detectorId = trackRef->DetectorId(); 
174           if (detectorId == 7)  refT0++;
175         }      
176       }
177     }
178   fTotalMult->Fill(numPrim);
179   fPrim->Fill(primaryC+primaryA);
180   fRef->Fill(refT0);
181   //     printf (" tracks %i primaries %d \n",nTracks, numPrim);
182  
183
184   // ESD information  
185   AliVEvent* event = InputEvent();
186   if (!event) {
187     Printf("ERROR: Could not retrieve event");
188     return;
189   }
190     
191   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
192
193   Bool_t eventTriggered = false;
194   ULong64_t triggerMask = esd->GetTriggerMask();
195   // definitions from p-p.cfg
196   ULong64_t spdFO = (1 << 14);
197   ULong64_t v0left = (1 << 11);
198   ULong64_t v0right = (1 << 12);
199
200   if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
201     eventTriggered == true;
202  
203     //= AliPWG0Helper::IsEventTriggered(esd, AliPWG0Helper::kMB1);
204   //  printf("!!!!! eventTriggered %i",eventTriggered);
205   // if(!eventTriggered) return;
206    Double_t besttimeC=99999.;
207   Double_t besttimeA=99999.;
208   
209 //  Float_t coefA = 0.891466;
210 //  Float_t coefC = 0.922772;
211   Float_t coefA = 1;
212   Float_t coefC = 1;
213   Float_t sumampA=0;
214   Float_t sumampC=0;
215   Int_t highA=0, highC=0;
216   const Double_t *amp = esd->GetT0amplitude();
217   const Double_t *time = esd->GetT0time();
218   Float_t vertex = esd->GetT0zVertex();
219   if(vertex<999) fVertex->Fill(vertex);
220   Float_t   mean = esd->GetT0();
221   fMean->Fill(mean);
222   
223   for (Int_t i=0; i<12; i++) {
224     sumampC += amp[i];  
225     fTime->Fill(amp[i],time[i]); 
226     fAmp->Fill(amp[i]); 
227     if(time[i]<besttimeC && time[i]>0) besttimeC=time[i]; //timeC
228   }
229   fOrC->Fill(besttimeC);
230   
231   for (Int_t i=12; i<24; i++){
232     sumampA += amp[i];  
233     fAmp->Fill(amp[i]); 
234     fTime->Fill(time[i]);
235     if(time[i]<besttimeA && time[i]>0) besttimeA=time[i]; //timeC
236   }
237
238   fOrA->Fill(besttimeA);
239   
240   fMultRecRealA ->Fill(primaryA, sumampA*coefA);
241   fMultRecRealC ->Fill(primaryC, sumampC*coefC);
242   fMultRecTotal->Fill(numPrim, sumampA*coefA +  sumampC*coefC);
243   fRec->Fill(Int_t(sumampA+sumampC));
244  
245   // Post output data.
246   PostData(1, fListOfHistos);
247 }      
248
249
250 void AliT0MultiplicityTask::Terminate(Option_t *) 
251 {
252   // Draw result to the screen
253   // Called once at the end of the query 
254   //  printf(" AliT0MultiplicityTask::Terminate ");
255   
256   fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
257   if (!fListOfHistos) {
258     Printf("ERROR: fListOfHistos not available");
259     return;
260   }
261   //  printf(" before cast ");
262   fOrA = dynamic_cast<TH1F*>(fListOfHistos->At(0));
263   fOrC = dynamic_cast<TH1F*>(fListOfHistos->At(1));
264   fMean  = dynamic_cast<TH1F*>(fListOfHistos->At(2));  
265   fVertex  = dynamic_cast<TH1F*>(fListOfHistos->At(3));
266   fAmp = dynamic_cast<TH1F*>(fListOfHistos->At(4));
267   fTime = dynamic_cast<TH1F*>(fListOfHistos->At(5));
268   fTotalMult = dynamic_cast<TH1F*>(fListOfHistos->At(6));
269   fMultRecTotal = dynamic_cast<TH2F*>(fListOfHistos->At(7));
270   fMultRecRealA = dynamic_cast<TH2F*>(fListOfHistos->At(8));
271   fMultRecRealC= dynamic_cast<TH2F*>(fListOfHistos->At(9));
272   fPrim= dynamic_cast<TH1F*>(fListOfHistos->At(10));
273   fRef= dynamic_cast<TH1F*>(fListOfHistos->At(11));
274   fRec= dynamic_cast<TH1F*>(fListOfHistos->At(12));
275
276   TFile fc("MultHist.root\n", "RECREATE");
277   //  printf(" File MultHist.root recreated\n");
278   fOrC->Write();
279   //  printf("1\n");
280   fOrA->Write();
281   //  printf("2\n");
282   fMean->Write();
283   //  printf("3\n");
284   fVertex->Write();
285   //  printf("4\n");
286   fAmp->Write();
287   //  printf("5\n");
288   fTime->Write();
289  fTotalMult->Write();
290  //  printf("6\n");
291  fMultRecTotal->Write();
292  //  printf("7\n");
293  fMultRecRealA->Write();
294  // printf("8\n");
295   fMultRecRealC->Write();
296   // printf("9\n");
297   fPrim->Write();
298   fRef->Write();
299   fRec->Write();
300   
301   fc.Close();
302  printf(" fc.Close()\n");
303
304 }