AliAnalysisTaskTriggerStudy: Added code for trigger overlap (venn-like) and loop...
[u/mrichter/AliRoot.git] / PWG0 / multPbPb / AliAnalysisTaskTriggerStudy.cxx
1 // AliAnalysisTaskTriggerStudy
2
3 // Author: Michele Floris, CERN
4 // TODO:
5 // - Add chi2/cluster plot for primary, secondaries and fakes
6
7
8 #include "AliAnalysisTaskTriggerStudy.h"
9 #include "AliESDInputHandler.h"
10 #include "AliHistoListWrapper.h"
11 #include "AliAnalysisManager.h"
12 #include "AliMCEvent.h"
13 #include "AliStack.h"
14 #include "TH1I.h"
15 #include "TH3D.h"
16 #include "AliMCParticle.h"
17 #include "AliGenEventHeader.h"
18 #include "AliESDCentrality.h"
19
20 #include <iostream>
21 #include "AliTriggerAnalysis.h"
22 #include "AliMultiplicity.h"
23 #include "TFile.h"
24 #include "AliLog.h"
25
26 using namespace std;
27
28 ClassImp(AliAnalysisTaskTriggerStudy)
29
30 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy()
31 : AliAnalysisTaskSE("TaskTriggerStudy"),
32   fESD(0),fHistoList(0),fIsMC(0),fTriggerAnalysis(0),fHistoSuffix("")
33 {
34   // constructor
35
36   DefineOutput(1, AliHistoListWrapper::Class());
37
38 }
39 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy(const char * name)
40   : AliAnalysisTaskSE(name),
41     fESD(0),fHistoList(0),fIsMC(0),fTriggerAnalysis(0),fHistoSuffix("")
42 {
43   //
44   // Standard constructur which should be used
45   //
46
47   DefineOutput(1, AliHistoListWrapper::Class());
48
49 }
50
51 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy(const AliAnalysisTaskTriggerStudy& obj) : 
52   AliAnalysisTaskSE(obj) ,fESD (0), fIsMC(0), fTriggerAnalysis(0),fHistoSuffix("")
53 {
54   //copy ctor
55   fESD = obj.fESD ;
56   fHistoList = obj.fHistoList;
57   fTriggerAnalysis = obj.fTriggerAnalysis;
58   fHistoSuffix = obj.fHistoSuffix;
59 }
60
61 AliAnalysisTaskTriggerStudy::~AliAnalysisTaskTriggerStudy(){
62   // destructor
63
64   if(!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
65     if(fHistoList) {
66       delete fHistoList;
67       fHistoList = 0;
68     }
69     if(fTriggerAnalysis) {
70       delete fTriggerAnalysis;
71       fHistoList = 0;
72     }
73   }
74   // Histo list should not be destroyed: fListWrapper is owner!
75
76 }
77 void AliAnalysisTaskTriggerStudy::UserCreateOutputObjects()
78 {
79   // Called once
80   fHistoList = new AliHistoListWrapper("histoList","histogram list for trigger studies");
81   fTriggerAnalysis = new AliTriggerAnalysis();
82 }
83
84
85 void AliAnalysisTaskTriggerStudy::UserExec(Option_t *)
86 {
87   // User code
88
89   /* PostData(0) is taken care of by AliAnalysisTaskSE */
90   PostData(1,fHistoList);
91
92   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
93   if (strcmp(fESD->ClassName(),"AliESDEvent")) {
94     AliFatal("Not processing ESDs");
95   }
96
97   // FIXME: two options here: either we add a loop setting the name of
98   // the histos with the trigger class in them (there may be a global
99   // SetHistoNamePrefix) or I put a cut to select only collision
100   // classes
101   
102   // get the multiplicity object
103   const AliMultiplicity* mult = fESD->GetMultiplicity();
104   Int_t ntracklets = mult->GetNumberOfTracklets();
105
106   // Reset histo suffix and fill reference histograms without any suffix
107   fHistoSuffix = "";
108   GetHistoTracklets("all","All events")->Fill(ntracklets);
109
110   // Fast or in the outer layer  
111   Int_t nFastOrOnline  = fTriggerAnalysis->SPDFiredChips(fESD, 1, 0, 2); // offline
112   Int_t nFastOrOffline = fTriggerAnalysis->SPDFiredChips(fESD, 0, 0, 2); // online
113   
114   Bool_t c0sm1 = nFastOrOffline >= 1;
115   Bool_t c0sm2 = nFastOrOffline >= 2;
116   Bool_t c0sm3 = nFastOrOffline >= 3;
117   Bool_t c0sm4 = nFastOrOffline >= 4;
118   Bool_t c0sm5 = nFastOrOffline >= 5;
119   
120   // V0 triggers
121   Bool_t c0v0A       = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0A);
122   Bool_t c0v0C       = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0C);
123   Bool_t v0AHW     = (fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
124   Bool_t v0CHW     = (fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
125
126   // TOF triggers 
127   // FIXME: implement real triggers
128   Bool_t c0OM2 = kFALSE;
129   Bool_t c0OM3 = kFALSE;
130
131   // Some macros for the online triggers
132   Bool_t cMBS2A = c0sm2 && c0v0A;
133   Bool_t cMBS2C = c0sm2 && c0v0C;
134   Bool_t cMBAC  = c0v0A && c0v0C;
135
136   FillTriggerOverlaps("All", "All Events", nFastOrOffline,c0v0A,c0v0C,c0OM2,c0OM3,cMBS2A,cMBS2C,cMBAC);
137   
138
139   // loop over trigger classes in the event
140   TObjArray * tokens = 0;
141   if(fIsMC) {
142     // in case of montecarlo I override the trigger class
143     tokens = new TObjArray;
144     tokens->SetOwner();
145     //    tokens->Add(new TObjString("CINT1B-ABCE-NOPF-ALL")); 
146     tokens->Add(new TObjString("MC")); 
147   }
148   else {  
149     TString trgClasses = fESD->GetFiredTriggerClasses();
150     tokens = trgClasses.Tokenize(" ");
151   }
152   TIter iter(tokens);
153     
154   while(TObjString * tok = (TObjString*) iter.Next()){
155     // clean up trigger name
156     TString trg = tok->GetString();
157     trg.Strip(TString::kTrailing, ' ');
158     trg.Strip(TString::kLeading, ' ');
159
160     fHistoSuffix = "_";
161     fHistoSuffix += trg;
162
163     // Fill histograms mismatchs
164     // TODO: check mismatch trigger class 
165     if(nFastOrOffline != nFastOrOnline) {
166       GetHistoTracklets("mismatchingFastOr", "Events where fast or offline differs from fast-or online")->Fill(ntracklets);
167     }
168     
169     if (c0v0A != v0AHW){
170       GetHistoTracklets("mismatchingV0A", "Events where V0A offline differs from V0A online")->Fill(ntracklets);
171     }
172     
173     if (c0v0C != v0CHW){
174       GetHistoTracklets("mismatchingV0C", "Events where V0C offline differs from V0C online")->Fill(ntracklets);
175     }    
176     
177     // Fill a tracklet histo for each trigger type
178     if(c0sm1)  GetHistoTracklets("c0sm1" ,"Events were trigger c0sm1 fired" )->Fill(ntracklets);
179     if(c0sm2)  GetHistoTracklets("c0sm2" ,"Events were trigger c0sm2 fired" )->Fill(ntracklets);
180     if(c0sm3)  GetHistoTracklets("c0sm3" ,"Events were trigger c0sm3 fired" )->Fill(ntracklets);
181     if(c0sm4)  GetHistoTracklets("c0sm4" ,"Events were trigger c0sm4 fired" )->Fill(ntracklets);
182     if(c0sm5)  GetHistoTracklets("c0sm5" ,"Events were trigger c0sm5 fired" )->Fill(ntracklets);
183     if(c0OM2)  GetHistoTracklets("c0OM2" ,"Events were trigger c0OM2 fired" )->Fill(ntracklets);
184     if(c0OM3)  GetHistoTracklets("c0OM3" ,"Events were trigger c0OM3 fired" )->Fill(ntracklets);
185     if(c0v0A)  GetHistoTracklets("c0v0A" ,"Events were trigger c0v0A fired" )->Fill(ntracklets);
186     if(c0v0C)  GetHistoTracklets("c0v0C" ,"Events were trigger c0v0C fired" )->Fill(ntracklets);
187     if(cMBS2A) GetHistoTracklets("cMBS2A","Events were trigger cMBS2A fired")->Fill(ntracklets);
188     if(cMBS2C) GetHistoTracklets("cMBS2C","Events were trigger cMBS2C fired")->Fill(ntracklets);
189     if(cMBAC ) GetHistoTracklets("cMBAC ","Events were trigger cMBAC  fired")->Fill(ntracklets);
190     //  if() GetHistoTracklets("","Events were trigger  fired");
191     delete tokens;
192   }
193     
194     // if (fIsMC) {
195     
196
197   //   if (!fMCEvent) {
198   //     AliError("No MC info found");
199   //   } else {
200       
201   //     //loop on the MC event
202   //     //      Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
203   //     Int_t offset    = fMCEvent->GetPrimaryOffset();
204   //     Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries()+offset;
205   //     for (Int_t ipart=offset; ipart<nMCTracks; ipart++) { 
206         
207   //    AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
208         
209   //    // We don't care about neutrals and non-physical primaries
210   //    if(mcPart->Charge() == 0) continue;
211
212   //    // FIXME: add kTransportBit (uncomment below)
213   //    if(!fMCEvent->Stack()->IsPhysicalPrimary(ipart)) continue;
214
215   //    //check if current particle is a physical primary
216   //    // Bool_t physprim=fMCEvent->IsPhysicalPrimary(label);
217   //    // if (!physprim) continue;
218   //    // if (!track) return kFALSE;
219   //    // Bool_t transported = mcPart->Particle()->TestBit(kTransportBit);
220   //    // if(!transported) return kFALSE;
221  
222   //    // Get MC vertex
223   //    //FIXME: which vertex do I take for MC?
224   //    TArrayF   vertex;
225   //    fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
226   //    Float_t zv = vertex[2];
227   //    //      Float_t zv = vtxESD->GetZ();
228   //    // Fill generated histo
229   //    hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv);
230         
231   //     }
232   //   }
233   // }
234   
235
236
237
238 }
239
240 void   AliAnalysisTaskTriggerStudy::Terminate(Option_t *){
241   // terminate
242   // Save output in a more friendly format
243   fHistoList = dynamic_cast<AliHistoListWrapper*> (GetOutputData(1));
244   if (!fHistoList){
245     Printf("ERROR: fHistoList not available");
246     return;
247   }
248   TFile * f = new TFile("trigger_study.root", "recreate");
249   fHistoList->GetList()->Write();
250   f->Close();
251
252 }
253
254 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoTracklets(const char * name, const char * title){
255   // Book histo of events vs ntracklets, if needed
256
257   TString hname = "hTracklets_";
258   hname+=name;  
259   hname+=fHistoSuffix;
260   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
261   
262   if(!h) {
263     AliInfo(Form("Booking histo %s",hname.Data()));
264     Bool_t oldStatus = TH1::AddDirectoryStatus();
265     TH1::AddDirectory(kFALSE);
266     h = new TH1F (hname.Data(), title, 50, 0.5, 200);
267     h->Sumw2();
268     h->SetXTitle("ntracklets");
269     fHistoList->GetList()->Add(h);
270     TH1::AddDirectory(oldStatus);
271   }
272   return h;
273 }
274
275 void AliAnalysisTaskTriggerStudy::FillTriggerOverlaps (const char * name, const char * title, 
276                                                        Int_t nFastOrOffline,
277                                                        Bool_t v0A, Bool_t v0C, Bool_t OM2, Bool_t OM3, 
278                                                        Bool_t cMBS2A,Bool_t cMBS2C, Bool_t cMBAC) {
279   //Fills a histo with the different trigger statistics in a venn like diagramm. Books it if needed.
280
281   // Get or book histo
282   TString hname = "hTrigStat_";
283   hname+=name;  
284   hname+=fHistoSuffix;
285   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
286   
287   if(!h) {
288     AliInfo(Form("Booking histo %s",hname.Data()));
289     Bool_t oldStatus = TH1::AddDirectoryStatus();
290     TH1::AddDirectory(kFALSE);
291     Int_t nbins = 14;
292     h = new TH1I (hname, title, nbins, 0.5, nbins+0.5);
293     fHistoList->GetList()->Add(h);
294     TH1::AddDirectory(oldStatus);
295   }
296
297   // we look at the combinations of 4 triggers
298   Int_t ibin = 1;
299   h->GetXaxis()->SetBinLabel(ibin++,"FO2 & V0A & V0C & OM2"); 
300   h->GetXaxis()->SetBinLabel(ibin++,"FO2 & V0C & OM2"); 
301   h->GetXaxis()->SetBinLabel(ibin++,"FO2 & V0A & OM2"); 
302   h->GetXaxis()->SetBinLabel(ibin++,"FO2 & V0A & V0C"); 
303   h->GetXaxis()->SetBinLabel(ibin++,"FO2 & V0A"); 
304   h->GetXaxis()->SetBinLabel(ibin++,"FO2 & V0C"); 
305   h->GetXaxis()->SetBinLabel(ibin++,"FO2 & OM2"); 
306   h->GetXaxis()->SetBinLabel(ibin++,"OM2 & V0A"); 
307   h->GetXaxis()->SetBinLabel(ibin++,"OM2 & V0C"); 
308   h->GetXaxis()->SetBinLabel(ibin++,"V0A & V0C"); 
309   h->GetXaxis()->SetBinLabel(ibin++,"FO2"); 
310   h->GetXaxis()->SetBinLabel(ibin++,"V0A"); 
311   h->GetXaxis()->SetBinLabel(ibin++,"V0C"); 
312   h->GetXaxis()->SetBinLabel(ibin++,"OM2"); 
313   
314   Bool_t fo2 = nFastOrOffline>=2;
315
316   if(fo2 && v0A && v0C && OM2)     {cout << "Bin5: " << ibin << endl;  h->Fill(1);}
317   if(fo2 && !v0A && v0C && OM2)    {cout << "Bin6: " << ibin << endl;  h->Fill(2);}
318   if(fo2 && v0A && !v0C && OM2)    {cout << "Bin7: " << ibin << endl;  h->Fill(3);}
319   if(fo2 && v0A && v0C && !OM2)    {cout << "Bin8: " << ibin << endl;  h->Fill(4);}
320   if(fo2 && v0A && !v0C && !OM2)   {cout << "Bin9: " << ibin << endl;  h->Fill(5);}
321   if(fo2 && !v0A && v0C && !OM2)   {cout << "Bin10: " << ibin << endl; h->Fill(6);}
322   if(fo2 && !v0A && !v0C && OM2)   {cout << "Bin11: " << ibin << endl; h->Fill(7);}
323   if(!fo2 && v0A && !v0C && OM2)   {cout << "Bin12: " << ibin << endl; h->Fill(8);}
324   if(!fo2 && !v0A && v0C && OM2)   {cout << "Bin13: " << ibin << endl; h->Fill(9);}
325   if(!fo2 && v0A && v0C && !OM2)   {cout << "Bin14: " << ibin << endl; h->Fill(10);}
326   if(fo2 && !v0A && !v0C && !OM2)  {cout << "Bin1: " << ibin << endl;  h->Fill(11);}
327   if(!fo2 && v0A && !v0C && !OM2)  {cout << "Bin2: " << ibin << endl;  h->Fill(12);}
328   if(!fo2 && !v0A && v0C && !OM2)  {cout << "Bin3: " << ibin << endl;  h->Fill(13);}
329   if(!fo2 && !v0A && !v0C && OM2)  {cout << "Bin4: " << ibin << endl;  h->Fill(14);}
330
331 }