Migrating PWG2/SPECTRA/Fit to new PWG structure
[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 "AliCentrality.h"
19
20 #include <iostream>
21 #include "AliTriggerAnalysis.h"
22 #include "AliMultiplicity.h"
23 #include "TFile.h"
24 #include "AliLog.h"
25 #include "AliESDtrackCuts.h"
26 #include "AliESDVZERO.h"
27 #include "TH2F.h"
28 #include "AliESDUtils.h"
29 #include "AliGenPythiaEventHeader.h"
30 #include "AliGenDPMjetEventHeader.h"
31
32 using namespace std;
33
34 ClassImp(AliAnalysisTaskTriggerStudy)
35
36 //const char * AliAnalysisTaskTriggerStudy::kVDNames[] = {"C0SM1","C0SM2","C0VBA","C0VBC","C0OM2"};       
37 const char * AliAnalysisTaskTriggerStudy::kVDNames[] = {"V0AND online","V0AND offline","Physics Selection", "Rec Candle", "Gen Canlde"};//,"C0OM2"};       
38
39 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy()
40 : AliAnalysisTaskSE("TaskTriggerStudy"),
41   fESD(0),fHistoList(0),fIsMC(0),fTriggerAnalysis(0),fHistoSuffix(""),fNTrackletsCut(1000000),fNTrackletsCutKine(100),fRejectBGWithV0(0)
42 {
43   // constructor
44
45   DefineOutput(1, AliHistoListWrapper::Class());
46
47 }
48 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy(const char * name)
49   : AliAnalysisTaskSE(name),
50     fESD(0),fHistoList(0),fIsMC(0),fTriggerAnalysis(0),fHistoSuffix(""),fNTrackletsCut(1000000),fNTrackletsCutKine(100),fRejectBGWithV0(0)
51 {
52   //
53   // Standard constructur which should be used
54   //
55
56   DefineOutput(1, AliHistoListWrapper::Class());
57
58 }
59
60 AliAnalysisTaskTriggerStudy::AliAnalysisTaskTriggerStudy(const AliAnalysisTaskTriggerStudy& obj) : 
61   AliAnalysisTaskSE(obj) ,fESD (0), fIsMC(0), fTriggerAnalysis(0),fHistoSuffix(""),fNTrackletsCut(1000000),fNTrackletsCutKine(100),fRejectBGWithV0(0)
62 {
63   //copy ctor
64   fESD = obj.fESD ;
65   fHistoList = obj.fHistoList;
66   fTriggerAnalysis = obj.fTriggerAnalysis;
67   fHistoSuffix = obj.fHistoSuffix;
68   fNTrackletsCut = obj.fNTrackletsCut;
69   fNTrackletsCutKine = obj.fNTrackletsCutKine;
70   fRejectBGWithV0 = obj.fRejectBGWithV0;
71 }
72
73 AliAnalysisTaskTriggerStudy::~AliAnalysisTaskTriggerStudy(){
74   // destructor
75
76   if(!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
77     if(fHistoList) {
78       delete fHistoList;
79       fHistoList = 0;
80     }
81     if(fTriggerAnalysis) {
82       delete fTriggerAnalysis;
83       fHistoList = 0;
84     }
85   }
86   // Histo list should not be destroyed: fListWrapper is owner!
87
88 }
89 void AliAnalysisTaskTriggerStudy::UserCreateOutputObjects()
90 {
91   // Called once
92   fHistoList = new AliHistoListWrapper("histoList","histogram list for trigger studies");
93   fTriggerAnalysis = new AliTriggerAnalysis();
94   if (fIsMC) fTriggerAnalysis->SetAnalyzeMC(1);
95 }
96
97
98 void AliAnalysisTaskTriggerStudy::UserExec(Option_t *)
99 {
100   // User code
101
102   // FIXME: make sure you have the right cuts here
103
104   /* PostData(0) is taken care of by AliAnalysisTaskSE */
105   PostData(1,fHistoList);
106
107   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
108   if (strcmp(fESD->ClassName(),"AliESDEvent")) {
109     AliFatal("Not processing ESDs");
110   }
111
112   // get the multiplicity object
113   const AliMultiplicity* mult = fESD->GetMultiplicity();
114   Int_t ntracklets = mult->GetNumberOfTracklets();
115   // Get Number of tracks
116   Int_t ntracks    = AliESDtrackCuts::GetReferenceMultiplicity(fESD,kTRUE); // tpc only
117   
118   // Get V0 Multiplicity
119   AliESDVZERO* esdV0 = fESD->GetVZEROData();
120   Float_t multV0A=esdV0->GetMTotV0A();
121   Float_t multV0C=esdV0->GetMTotV0C();
122   Float_t dummy = 0;  
123   Float_t multV0 =  AliESDUtils::GetCorrV0(fESD, dummy);
124
125   // Get number of clusters in layer 1
126   Float_t outerLayerSPD = mult->GetNumberOfITSClusters(1);  
127   Float_t innerLayerSPD = mult->GetNumberOfITSClusters(0);  
128   Float_t totalClusSPD = outerLayerSPD+innerLayerSPD;
129
130   //  if ( !fIsMC &&(!(fESD->IsTriggerClassFired("CMBS2A-B-NOPF-ALL")|| fESD->IsTriggerClassFired("CMBS2C-B-NOPF-ALL") || fESD->IsTriggerClassFired("CMBAC-B-NOPF-ALL")) )) return;
131   //  if ( !fIsMC &&(!(fESD->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")))) return;
132   if ( !fIsMC &&(!(fESD->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))) return;
133   GetHistoSPD1  ("All", "All events before any selection")->Fill(outerLayerSPD);
134   GetHistoV0M   ("All", "All events before any selection")->Fill(multV0);
135
136   // V0 triggers
137   Bool_t c0v0A       = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0A);
138   Bool_t c0v0C       = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0C);
139   // Bool_t v0AHW     = (fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
140   // Bool_t v0CHW     = (fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kTRUE) == AliTriggerAnalysis::kV0BB);// should replay hw trigger
141   Bool_t v0AHW = fTriggerAnalysis->EvaluateTrigger(fESD, AliTriggerAnalysis::kCTPV0A);
142   Bool_t v0CHW = fTriggerAnalysis->EvaluateTrigger(fESD, AliTriggerAnalysis::kCTPV0C);
143   Bool_t v0ABG = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0ABG);
144   Bool_t v0CBG = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0CBG);
145   Bool_t v0BG  = v0ABG || v0CBG;
146
147   // At least one track in eta < 0.8 with pt > 0.5
148   // MC Checks
149   Bool_t atLeast1Track = kFALSE;
150   Bool_t isSD = kFALSE;
151   Bool_t isND = kFALSE;
152
153   if(fIsMC) {
154     if (!fMCEvent) {
155       AliError("No MC info found");
156     } else {
157       Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
158       Int_t nPhysicalPrimaries = 0;
159       for (Int_t ipart=0; ipart<nMCTracks; ipart++) {   
160         AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
161         // We don't care about neutrals and non-physical primaries
162         if(mcPart->Charge() == 0) continue;
163
164         //check if current particle is a physical primary
165         if(!fMCEvent->IsPhysicalPrimary(ipart)) continue;
166         if(mcPart->Pt()<0.5) continue;
167         if(TMath::Abs(mcPart->Eta())>0.8) continue;
168         atLeast1Track = kTRUE;
169         break;
170       }
171       
172       AliGenPythiaEventHeader * headPy  = 0;
173       AliGenDPMjetEventHeader * headPho = 0;
174       AliGenEventHeader * htmp = fMCEvent->GenEventHeader();
175       if(!htmp) {
176         AliFatal("Cannot Get MC Header!!");
177         return;
178       }
179       if( TString(htmp->IsA()->GetName()) == "AliGenPythiaEventHeader") {
180         headPy =  (AliGenPythiaEventHeader*) htmp;
181       } else if (TString(htmp->IsA()->GetName()) == "AliGenDPMjetEventHeader") {
182         headPho = (AliGenDPMjetEventHeader*) htmp;
183       } else {
184         AliFatal("Unknown header");     
185       }
186       if(headPy)   {
187         //    cout << "Process: " << headPy->ProcessType() << endl;
188         if(headPy->ProcessType() == 92 || headPy->ProcessType() == 93) {
189           isSD = kTRUE; // is single difractive
190         }
191         if(headPy->ProcessType() != 92 && headPy->ProcessType() != 93 && headPy->ProcessType() != 94) {     
192           isND = kTRUE; // is non-diffractive
193         }
194         
195       } else if (headPho) {
196         if(headPho->ProcessType() == 5 || headPho->ProcessType() == 6 ) {
197           isSD = kTRUE;
198         }       
199         if(headPho->ProcessType() != 5 && headPho->ProcessType() != 6  && headPho->ProcessType() != 7 ) {
200           isND = kTRUE;
201         }       
202       }
203       
204     }
205   }
206   
207   static AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
208   // cuts->SetPtRange(0.5,10000);
209   // cuts->SetEtaRange(-0.8, 0.8);
210   Int_t ntracksLoop = fESD->GetNumberOfTracks();
211
212   Bool_t atLeast1TrackESD = kFALSE;
213   for (Int_t iTrack = 0; iTrack<ntracksLoop; iTrack++) {    
214     AliESDtrack * track = dynamic_cast<AliESDtrack*>(fESD->GetTrack(iTrack));
215     // for each track
216     
217     // track quality cuts
218     if(!cuts->AcceptTrack(track)) continue;
219     if(track->Pt()<0.5) continue;
220     if(TMath::Abs(track->Eta())>0.8) continue;
221
222     atLeast1TrackESD = kTRUE;
223     break;
224   }
225
226
227
228   // Physics selection
229   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
230
231   Bool_t vdArray[kNVDEntries];
232   vdArray[kVDV0ANDOnline]      = v0AHW && v0CHW && !v0BG ;
233   vdArray[kVDV0ANDOffline]     = c0v0A && c0v0C && isSelected;
234   vdArray[kVDRecCandle]        = atLeast1TrackESD;
235   vdArray[kVDPhysSel]          = isSelected;
236   vdArray[kVDGenCandle]        = atLeast1Track;
237
238   FillTriggerOverlaps("All", "All Events",vdArray);
239   if(!isSD)   FillTriggerOverlaps("NSD", "NSD Events",vdArray);
240
241
242   if(!isSelected) return;
243
244
245   GetHistoSPD1  ("AllPSNoPrino", "All events after physsel (no ZDC time)")->Fill(outerLayerSPD);
246   GetHistoV0M   ("AllPSNoPrino", "All events after physsel (no ZDC time)")->Fill(multV0);
247
248
249   // Francesco's cuts
250   // const AliESDVertex * vtxESDTPC= fESD->GetPrimaryVertexTPC(); 
251   // if(vtxESDTPC->GetNContributors()<1) return;
252   // if (vtxESDTPC->GetNContributors()<(-10.+0.25*fESD->GetMultiplicity()->GetNumberOfITSClusters(0)))     return;
253   // const AliESDVertex * vtxESDSPD= fESD->GetPrimaryVertexSPD(); 
254   // Float_t tpcContr=vtxESDTPC->GetNContributors();
255
256   
257
258   // GetT0 Stuff
259   const Double32_t *meanT0   = fESD->GetT0TOF();
260   const Double32_t  meanT0A  = 0.001* meanT0[1];
261   const Double32_t  meanT0C  = 0.001* meanT0[2];
262   const Double32_t  meanT0AC = 0.001* meanT0[0];
263   Double32_t  T0Vertex = fESD->GetT0zVertex();
264   //  Double32_t  *ampT0 =Esdevent ->GetT0amplitude();
265
266 //   cut1yesF = ( (meanC < 95. && meanA < 95.) && (meanC < -2.) ) && francescoscut
267 // cut1notF = ( (meanC < 95. && meanA < 95.) && (meanC < -2.) ) && ! francescoscut
268 // cut2 = ( (meanC < 95. && meanA < 95.) && ( (meanC-meanA) <=-0.7) && meanC > -2) )
269 // cut3 = ( (meanC < 95. && meanA < 95.) && ( (meanC-meanA) < 0.7 && (meanC-meanA) > -0.7 ) )
270 //   cut4 = ( (meanC < 95. && meanA < 95.) && (meanA < -2.)
271
272   Bool_t cut1T0 =  ( (meanT0C < 95. && meanT0A < 95.) && (meanT0C < -2.) );
273   Bool_t cut2T0 = ( (meanT0C < 95. && meanT0A < 95.) && ( (meanT0C-meanT0A) <=-0.7) && meanT0C > -2) ;
274
275   if(ntracklets > fNTrackletsCut) return;
276
277   // Reset histo suffix and fill reference histograms without any suffix
278   fHistoSuffix = "";
279
280   // Fast or in the outer layer  
281   Int_t nFastOrOnline  = fTriggerAnalysis->SPDFiredChips(fESD, 1, 0, 2); // offline
282   Int_t nFastOrOffline = fTriggerAnalysis->SPDFiredChips(fESD, 0, 0, 2); // online
283   
284   Bool_t c0sm1 = nFastOrOffline >= 1;
285   Bool_t c0sm2 = nFastOrOffline >= 2;
286   Bool_t c0sm3 = nFastOrOffline >= 3;
287   Bool_t c0sm4 = nFastOrOffline >= 4;
288   Bool_t c0sm5 = nFastOrOffline >= 5;
289   
290
291   // TOF triggers 
292   // FIXME: move to triggeranalysis?
293   AliESDHeader*h = fESD->GetHeader(); // taken the header from AliESDEvent 
294   Bool_t c0OM2 = h->IsTriggerInputFired("0OM2"); // thr >= 2 (input 19)
295   Bool_t c0OM3 = h->IsTriggerInputFired("0OM3"); // thr >= 3 (input 20)
296
297   // ZDC triggers
298   Bool_t zdcA    = kFALSE;
299   Bool_t zdcC    = kFALSE;
300   Bool_t zdcBar  = kFALSE;
301   Bool_t zdcTime = fTriggerAnalysis->ZDCTimeTrigger(fESD);
302
303   if (!fIsMC) {
304     // If it's data, we use the TDCs
305     zdcA   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kASide, kTRUE, kFALSE) ; 
306     zdcC   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kCSide, kTRUE, kFALSE) ;                 
307     zdcBar = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kCentralBarrel) ;                                
308   } else {
309     // If it's MC, we use the energy
310     Double_t minEnergy = 0;
311     AliESDZDC *esdZDC = fESD->GetESDZDC();    
312     Double_t  zNCEnergy = esdZDC->GetZDCN1Energy();
313     Double_t  zPCEnergy = esdZDC->GetZDCP1Energy();
314     Double_t  zNAEnergy = esdZDC->GetZDCN2Energy();
315     Double_t  zPAEnergy = esdZDC->GetZDCP2Energy();
316     // zdcA = (zNAEnergy>minEnergy || zPAEnergy>minEnergy);
317     // zdcC = (zNCEnergy>minEnergy || zPCEnergy>minEnergy);
318     zdcA = (zNAEnergy>minEnergy);
319     zdcC = (zNCEnergy>minEnergy);
320   }
321
322
323   // Some macros for the online triggers
324   Bool_t cMBS2A = c0sm2 && c0v0A;
325   Bool_t cMBS2C = c0sm2 && c0v0C;
326   Bool_t cMBAC  = c0v0A && c0v0C;
327   
328
329
330
331
332   // Plots requested by Jurgen on 18/11/2010 + later updates (including plots for the note)
333   // FIXME: will skip everything else
334
335
336   if(zdcTime) {
337     GetHistoSPD1  ("ZDCTIME", "ZDC Time Cut")->Fill(outerLayerSPD);
338     GetHistoTracks("ZDCTIME", "ZDC Time Cut")->Fill(ntracks);
339     GetHistoV0M   ("ZDCTIME", "ZDC Time Cut")->Fill(multV0);
340   }
341
342   if(zdcTime && ntracks > 1) {
343     GetHistoSPD1  ("ZDCTIME1TRACK", "ZDC Time Cut & 1 Track")->Fill(outerLayerSPD);
344     GetHistoTracks("ZDCTIME1TRACK", "ZDC Time Cut & 1 Track")->Fill(ntracks);
345     GetHistoV0M   ("ZDCTIME1TRACK", "ZDC Time Cut & 1 Track")->Fill(multV0);
346   }
347
348   // GetHistoSPD1  ("PhysSel", "All events after physics selection and Francesco's cut")->Fill(outerLayerSPD);
349   // GetHistoTracks("PhysSel", "All events after physics selection and Francesco's cut")->Fill(ntracks);
350   // GetHistoV0M   ("PhysSel", "All events after physics selection and Francesco's cut")->Fill(multV0);
351   if(c0v0A && c0v0C) {
352     GetHistoSPD1  ("V0AND", "V0A & V0C")->Fill(outerLayerSPD);
353     GetHistoTracks("V0AND", "V0A & V0C")->Fill(ntracks);
354     GetHistoV0M   ("V0AND", "V0A & V0C")->Fill(multV0);
355   }
356   if((c0v0A && !c0v0C) || (!c0v0A && c0v0C)) {
357     GetHistoSPD1  ("V0ONLYONE", "(V0A & !V0C) || (!V0A & V0C)")->Fill(outerLayerSPD);
358     GetHistoTracks("V0ONLYONE", "(V0A & !V0C) || (!V0A & V0C)")->Fill(ntracks);
359     GetHistoV0M   ("V0ONLYONE", "(V0A & !V0C) || (!V0A & V0C)")->Fill(multV0);
360   }
361   if(zdcA && zdcC) {
362     GetHistoSPD1  ("ZDCAND", "ZDCA & ZDCC")->Fill(outerLayerSPD);
363     GetHistoTracks("ZDCAND", "ZDCA & ZDCC")->Fill(ntracks);
364     GetHistoV0M   ("ZDCAND", "ZDCA & ZDCC")->Fill(multV0);
365   }
366   if((c0v0A && c0v0C) && !(zdcA && zdcC)) {
367     GetHistoSPD1  ("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(outerLayerSPD);
368     GetHistoTracks("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(ntracks);
369     GetHistoV0M   ("V0ANDNOTZDCAND", "(V0A & V0C) & !(ZDCA & ZDCC)")->Fill(multV0);
370   }
371   if((c0v0A && c0v0C) && !zdcA && !zdcC) {
372     GetHistoSPD1  ("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(outerLayerSPD);
373     GetHistoTracks("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(ntracks);
374     GetHistoV0M   ("V0ANDNOTZDC", "(V0A & V0C) & !ZDCA & !ZDCC)")->Fill(multV0);
375 }
376   if((c0v0A && c0v0C) && ((!zdcA && zdcC) || (zdcA && !zdcC))) {
377     GetHistoSPD1  ("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(outerLayerSPD);
378     GetHistoTracks("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(ntracks);
379     GetHistoV0M   ("V0ANDONEZDC", "(V0A & V0C) &  ((!ZDCA && ZDCC) || (ZDCA && !ZDCC)))")->Fill(multV0);
380   }
381
382   if(((c0v0A && !c0v0C) || (!c0v0A && c0v0C)) && (zdcA && zdcC)) {
383     GetHistoSPD1  ("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(outerLayerSPD);
384     GetHistoTracks("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(ntracks);
385     GetHistoV0M   ("V0ONEZDCBOTH", "((V0A && !V0C) || (!V0A && V0C)) && (ZDCA && ZDCC)")->Fill(multV0);
386   }
387
388   if((c0v0A && c0v0C) && (zdcA && zdcC)) {
389     GetHistoSPD1  ("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(outerLayerSPD);
390     GetHistoTracks("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(ntracks);
391     GetHistoV0M   ("V0ANDZDCAND", "(V0A & V0C) & (ZDCA & ZDCC)")->Fill(multV0);
392   }
393
394   if((c0v0A && zdcA && !zdcC && !c0v0C) || (c0v0C && zdcC && !zdcA && !c0v0A)) {
395     GetHistoSPD1  ("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(outerLayerSPD);
396     GetHistoTracks("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(ntracks);
397     GetHistoV0M   ("OneSided", "(V0A & ZDCA & !ZDCC & !V0C) || (V0C & ZDCC & !ZDCA & !V0A)")->Fill(multV0);
398   }
399
400   // GetHistoCorrelationSPDTPCVz("All", "After physics selection and Francesco's cut")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
401   // if(cut1T0)   GetHistoCorrelationSPDTPCVz("Cut1T0", "T0 Cut 1")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
402   // if(cut2T0)   GetHistoCorrelationSPDTPCVz("Cut2T0", "T0 Cut 2")->Fill(vtxESDSPD->GetZ(),vtxESDTPC->GetZ());
403
404   // GetHistoCorrelationContrTPCSPDCls("All", "After physics selection and Francesco's cut")->Fill(totalClusSPD,tpcContr);
405   // if(cut1T0)   GetHistoCorrelationContrTPCSPDCls("Cut1T0", "T0 Cut 1")->Fill(totalClusSPD,tpcContr);
406   // if(cut2T0)   GetHistoCorrelationContrTPCSPDCls("Cut2T0", "T0 Cut 2")->Fill(totalClusSPD,tpcContr);
407
408   // GetHistoCorrelationTrackletsSPDCls("All", "After physics selection and Francesco's cut")->Fill(totalClusSPD,ntracklets);
409   // if(cut1T0)   GetHistoCorrelationTrackletsSPDCls("Cut1T0", "T0 Cut 1")->Fill(totalClusSPD,ntracklets);
410   // if(cut2T0)   GetHistoCorrelationTrackletsSPDCls("Cut2T0", "T0 Cut 2")->Fill(totalClusSPD,ntracklets);
411
412
413   return; // FIXME
414
415
416   // Reject background
417   if (v0BG && fRejectBGWithV0) {
418     cout << "Rejection BG" << endl;
419     
420     return;
421   }
422   // Fill global histos
423   GetHistoTracklets("all","All events")->Fill(ntracklets);
424   FillTriggerOverlaps("All", "All Events",vdArray);
425   
426   // Fill some combination of trigger classes
427   Bool_t cmbs1aOnline = fESD->IsTriggerClassFired("CMBS1A-B-NOPF-ALL");
428   Bool_t cmbs1cOnline = fESD->IsTriggerClassFired("CMBS1C-B-NOPF-ALL");
429   Bool_t cmbacOnline  = fESD->IsTriggerClassFired("CMBAC-B-NOPF-ALL");
430   
431   Bool_t twoOutOfThree = kFALSE;
432   if (cmbs1aOnline || cmbs1cOnline ||cmbacOnline) twoOutOfThree = kTRUE;
433   if (twoOutOfThree) GetHistoTracklets("All_TwoOutOfThree" ,"Events 2-out-of-3 online" )->Fill(ntracklets);
434
435
436   // loop over trigger classes in the event
437   TObjArray * tokens = 0;
438   if(fIsMC) {
439     // in case of montecarlo I override the trigger class
440     tokens = new TObjArray;
441     tokens->SetOwner();
442     //    tokens->Add(new TObjString("CINT1B-ABCE-NOPF-ALL")); 
443     tokens->Add(new TObjString("MC")); 
444   }
445   else {  
446     TString trgClasses = fESD->GetFiredTriggerClasses();
447     tokens = trgClasses.Tokenize("  ");
448   }
449   TIterator  * iter = (TIterator*) tokens->MakeIterator();
450   
451   TString classes = fESD->GetFiredTriggerClasses();
452
453   // if (classes.Contains("SMH")) {
454   //    tokens->Print();
455   //    cout << classes.Data() << endl;
456   // }
457   //  iter->Reset();
458   //Int_t itoken = 0;
459   TObjString * tok=0;
460   while((tok = (TObjString*) iter->Next())){
461     // clean up trigger name
462     TString trg = tok->GetString();
463     trg.Strip(TString::kTrailing, ' ');
464     trg.Strip(TString::kLeading, ' ');
465
466     fHistoSuffix = "_";
467     fHistoSuffix += trg;
468
469     //    cout << itoken++ << " " << trg.Data() << endl;
470     // continue;
471     //    if (trg.Contains("SMH")) cout << itoken++ << " " << trg.Data() << endl;
472     
473     // Fill tracklets 
474     GetHistoTracklets("All" ,"Events no offline trigger" )->Fill(ntracklets);    
475
476
477     // Fill histograms mismatchs
478     // TODO: check mismatch trigger class 
479     if(nFastOrOffline != nFastOrOnline) {
480       GetHistoTracklets("mismatchingFastOr", "Events where fast or offline differs from fast-or online")->Fill(ntracklets);
481     }
482     
483     if (c0v0A != v0AHW){
484       GetHistoTracklets("mismatchingV0A", "Events where V0A offline differs from V0A online")->Fill(ntracklets);
485     }
486     
487     if (c0v0C != v0CHW){
488       GetHistoTracklets("mismatchingV0C", "Events where V0C offline differs from V0C online")->Fill(ntracklets);
489     }    
490     
491     // Fill a tracklet histo for each trigger type
492     if(c0sm1)  GetHistoTracklets("c0sm1" ,"Events were trigger c0sm1 fired" )->Fill(ntracklets);
493     if(c0sm2)  GetHistoTracklets("c0sm2" ,"Events were trigger c0sm2 fired" )->Fill(ntracklets);
494     if(c0sm3)  GetHistoTracklets("c0sm3" ,"Events were trigger c0sm3 fired" )->Fill(ntracklets);
495     if(c0sm4)  GetHistoTracklets("c0sm4" ,"Events were trigger c0sm4 fired" )->Fill(ntracklets);
496     if(c0sm5)  GetHistoTracklets("c0sm5" ,"Events were trigger c0sm5 fired" )->Fill(ntracklets);
497     if(c0OM2)  GetHistoTracklets("c0OM2" ,"Events were trigger c0OM2 fired" )->Fill(ntracklets);
498     if(c0OM3)  GetHistoTracklets("c0OM3" ,"Events were trigger c0OM3 fired" )->Fill(ntracklets);
499     if(c0v0A)  GetHistoTracklets("c0v0A" ,"Events were trigger c0v0A fired" )->Fill(ntracklets);
500     if(c0v0C)  GetHistoTracklets("c0v0C" ,"Events were trigger c0v0C fired" )->Fill(ntracklets);
501     if(cMBS2A) GetHistoTracklets("cMBS2A","Events were trigger cMBS2A fired")->Fill(ntracklets);
502     if(cMBS2C) GetHistoTracklets("cMBS2C","Events were trigger cMBS2C fired")->Fill(ntracklets);
503     if(cMBAC ) GetHistoTracklets("cMBAC", "Events were trigger cMBAC  fired")->Fill(ntracklets);
504     if(zdcA )  GetHistoTracklets("cZDCA", "Events were trigger cZDCA  fired")->Fill(ntracklets);
505     if(zdcC )  GetHistoTracklets("cZDCC", "Events were trigger cZDCC  fired")->Fill(ntracklets);    
506     if(!zdcA && !zdcC )  GetHistoTracklets("NoZDC", "Events were zdc trigger don't  fired")->Fill(ntracklets);
507     if( (zdcA && !zdcC) || (!zdcA && zdcC) )  GetHistoTracklets("OneSideZDC", "Events were only 1 ZDC  fired")->Fill(ntracklets);
508     if( (zdcA && zdcC) )  GetHistoTracklets("TwoSideZDC", "Events were both  ZDC  fired")->Fill(ntracklets);
509
510     if(twoOutOfThree) {
511       if(!zdcA && !zdcC )  GetHistoTracklets("NoZDC_TwoOutOfThree", "Events were zdc trigger don't  fired")->Fill(ntracklets);
512       if( (zdcA && !zdcC) || (!zdcA && zdcC) )  GetHistoTracklets("OneSideZDC_TwoOutOfThree", "Events were only 1 ZDC  fired")->Fill(ntracklets);
513       if( (zdcA && zdcC) )  GetHistoTracklets("TwoSideZDC_TwoOutOfThree", "Events were both  ZDC  fired")->Fill(ntracklets);
514     }
515     if(zdcBar) GetHistoTracklets("cZDCBar","Events were trigger cZDCB  fired")->Fill(ntracklets);
516     //  if() GetHistoTracklets("","Events were trigger  fired");
517     
518     // Fill trigger overlaps
519     FillTriggerOverlaps("All", "All Events in trigger class",vdArray);
520
521
522     // Fill kinematic variables for peripheral events
523     static AliESDtrackCuts * cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(); // FIXME: change this ?
524     if (ntracklets < fNTrackletsCutKine) {
525       // 1. Loop on tracklets
526       for(Int_t itracklet = 0; itracklet < ntracklets; itracklet++){      
527         GetHistoEta("All", Form("Tracklet #eta distribution, ntracklets < %d",fNTrackletsCutKine))->Fill(mult->GetEta(itracklet));
528       }
529       // 2. loop on tracks
530       TObjArray * goodTracks = cuts->GetAcceptedTracks(fESD);
531       TIterator * iterTracks = goodTracks->MakeIterator();
532       AliESDtrack * esdTrack = 0;
533       while ((esdTrack = (AliESDtrack*) iterTracks->Next())) {
534         GetHistoPt("All", Form("Tracklet p_{T} distribution, ntracklets < %d",fNTrackletsCutKine))->Fill(esdTrack->Pt());       
535       }
536       delete goodTracks;
537     }
538     
539
540   }
541   delete tokens;
542     
543     // if (fIsMC) {
544     
545
546   //   if (!fMCEvent) {
547   //     AliError("No MC info found");
548   //   } else {
549       
550   //     //loop on the MC event
551   //     //      Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
552   //     Int_t offset    = fMCEvent->GetPrimaryOffset();
553   //     Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries()+offset;
554   //     for (Int_t ipart=offset; ipart<nMCTracks; ipart++) { 
555         
556   //    AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
557         
558   //    // We don't care about neutrals and non-physical primaries
559   //    if(mcPart->Charge() == 0) continue;
560
561   // PHYSICAL PRIMARY
562   //    // Get MC vertex
563     //  TArrayF   vertex;
564   //    fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
565   //    Float_t zv = vertex[2];
566   //    //      Float_t zv = vtxESD->GetZ();
567   //    // Fill generated histo
568   //    hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv);
569         
570   //     }
571   //   }
572   // }
573   
574
575
576
577 }
578
579 void   AliAnalysisTaskTriggerStudy::Terminate(Option_t *){
580   // terminate
581   // Save output in a more friendly format
582   fHistoList = dynamic_cast<AliHistoListWrapper*> (GetOutputData(1));
583   if (!fHistoList){
584     Printf("ERROR: fHistoList not available");
585     return;
586   }
587   TFile * f = new TFile("trigger_study.root", "recreate");
588   fHistoList->GetList()->Write();
589   f->Close();
590
591 }
592
593 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoTracklets(const char * name, const char * title){
594   // Book histo of events vs ntracklets, if needed
595
596   TString hname = "hTracklets_";
597   hname+=name;  
598   hname+=fHistoSuffix;
599   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
600   
601   if(!h) {
602     AliInfo(Form("Booking histo %s",hname.Data()));
603     Bool_t oldStatus = TH1::AddDirectoryStatus();
604     TH1::AddDirectory(kFALSE);
605     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
606     h = new TH1F (hname.Data(), title, 4000, -0.5, 4000-0.5);
607     h->Sumw2();
608     h->SetXTitle("ntracklets");
609     fHistoList->GetList()->Add(h);
610     TH1::AddDirectory(oldStatus);
611   }
612   return h;
613 }
614
615
616 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoSPD1(const char * name, const char * title){
617   // Book histo of events vs ntracklets, if needed
618
619   TString hname = "hSPD1_";
620   hname+=name;  
621   hname+=fHistoSuffix;
622   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
623   
624   // const Int_t nbin = 118;
625   // Float_t bins[nbin] = {0};
626   // bins[0] = 0;
627   // for(Int_t ibin = 1; ibin <= nbin; ibin++){
628   //   if (ibin < 100) 
629   //     bins[ibin] = bins [ibin-1]+1; 
630   //   else if (ibin < 1000) 
631   //     bins[ibin] = bins [ibin-1]+100; 
632   //   else if (ibin < 10000) 
633   //     bins[ibin] = bins [ibin-1]+1000; 
634   // }
635   
636
637   if(!h) {
638     AliInfo(Form("Booking histo %s",hname.Data()));
639     Bool_t oldStatus = TH1::AddDirectoryStatus();
640     TH1::AddDirectory(kFALSE);
641     //    h = new TH1F (hname.Data(), title, nbin, bins);
642     //h = new TH1F (hname.Data(), title, 100, -0.5, 100-0.5);
643     h = new TH1F (hname.Data(), title, 10000, -0.5, 10000-0.5);
644     h->Sumw2();
645     h->SetXTitle("SPD1");
646     fHistoList->GetList()->Add(h);
647     TH1::AddDirectory(oldStatus);
648   }
649   return h;
650 }
651
652 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoV0M(const char * name, const char * title){
653   // Book histo of events vs ntracklets, if needed
654
655   TString hname = "hV0M_";
656   hname+=name;  
657   hname+=fHistoSuffix;
658   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
659   
660   if(!h) {
661     AliInfo(Form("Booking histo %s",hname.Data()));
662     Bool_t oldStatus = TH1::AddDirectoryStatus();
663     TH1::AddDirectory(kFALSE);
664     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
665     h = new TH1F (hname.Data(), title, 500, -0.5, 500-0.5);
666     h->Sumw2();
667     h->SetXTitle("V0M");
668     fHistoList->GetList()->Add(h);
669     TH1::AddDirectory(oldStatus);
670   }
671   return h;
672 }
673
674
675 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoTracks(const char * name, const char * title){
676   // Book histo of events vs ntracklets, if needed
677
678   TString hname = "hTPCTracks_";
679   hname+=name;  
680   hname+=fHistoSuffix;
681   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
682   
683   if(!h) {
684     AliInfo(Form("Booking histo %s",hname.Data()));
685     Bool_t oldStatus = TH1::AddDirectoryStatus();
686     TH1::AddDirectory(kFALSE);
687     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
688     h = new TH1F (hname.Data(), title, 50, -0.5, 50-0.5);
689     h->Sumw2();
690     h->SetXTitle("ntracks");
691     fHistoList->GetList()->Add(h);
692     TH1::AddDirectory(oldStatus);
693   }
694   return h;
695 }
696
697
698 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoEta(const char * name, const char * title){
699   // Book histo of events vs ntracklets, if needed
700
701   TString hname = "hEta_";
702   hname+=name;  
703   hname+=fHistoSuffix;
704   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
705   
706   if(!h) {
707     AliInfo(Form("Booking histo %s",hname.Data()));
708     Bool_t oldStatus = TH1::AddDirectoryStatus();
709     TH1::AddDirectory(kFALSE);
710     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
711     h = new TH1F (hname.Data(), title, 20, -2, 2);
712     h->Sumw2();
713     h->SetXTitle("#eta");
714     fHistoList->GetList()->Add(h);
715     TH1::AddDirectory(oldStatus);
716   }
717   return h;
718 }
719
720 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoPt(const char * name, const char * title){
721   // Book histo of pt distribution, if needed
722
723   TString hname = "hPt_";
724   hname+=name;  
725   hname+=fHistoSuffix;
726   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
727   
728   if(!h) {
729     AliInfo(Form("Booking histo %s",hname.Data()));
730     Bool_t oldStatus = TH1::AddDirectoryStatus();
731     TH1::AddDirectory(kFALSE);
732     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
733     h = new TH1F (hname.Data(), title, 100, -0.5, 1-0.5);
734     h->Sumw2();
735     h->SetXTitle("p_{T}");
736     fHistoList->GetList()->Add(h);
737     TH1::AddDirectory(oldStatus);
738   }
739   return h;
740 }
741
742 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationSPDTPCVz(const char * name, const char * title){
743   // Book histo of events vs ntracklets, if needed
744
745   TString hname = "hTPCvsSPD_";
746   hname+=name;  
747   hname+=fHistoSuffix;
748   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
749   
750   if(!h) {
751     AliInfo(Form("Booking histo %s",hname.Data()));
752     Bool_t oldStatus = TH1::AddDirectoryStatus();
753     TH1::AddDirectory(kFALSE);
754     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
755     h = new TH2F (hname.Data(), title, 80, -20, 20,  80, -20, 20);
756     h->Sumw2();
757     h->SetXTitle("SPD Vz");
758     h->SetYTitle("TPC Vz");
759     fHistoList->GetList()->Add(h);
760     TH1::AddDirectory(oldStatus);
761   }
762   return h;
763 }
764
765 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationContrTPCSPDCls(const char * name, const char * title){
766   // Book histo of events vs ntracklets, if needed
767
768   TString hname = "hContrTPCvsSPDCls_";
769   hname+=name;  
770   hname+=fHistoSuffix;
771   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
772   
773   if(!h) {
774     AliInfo(Form("Booking histo %s",hname.Data()));
775     Bool_t oldStatus = TH1::AddDirectoryStatus();
776     TH1::AddDirectory(kFALSE);
777     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
778     h = new TH2F (hname.Data(), title, 1000, 0, 9000,  1000, 0, 5000);
779     h->Sumw2();
780     h->SetXTitle("SPD Clusters");
781     h->SetYTitle("TPC V Contributors");
782     fHistoList->GetList()->Add(h);
783     TH1::AddDirectory(oldStatus);
784   }
785   return h;
786 }
787 TH1 *   AliAnalysisTaskTriggerStudy::GetHistoCorrelationTrackletsSPDCls(const char * name, const char * title){
788   // Book histo of events vs ntracklets, if needed
789
790   TString hname = "hTrackletsvsSPDCls_";
791   hname+=name;  
792   hname+=fHistoSuffix;
793   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
794   
795   if(!h) {
796     AliInfo(Form("Booking histo %s",hname.Data()));
797     Bool_t oldStatus = TH1::AddDirectoryStatus();
798     TH1::AddDirectory(kFALSE);
799     // h = new TH1F (hname.Data(), title, 1000, -0.5, 10000-0.5);
800     h = new TH2F (hname.Data(), title, 1000, 0, 9000,  1000, 0, 5000);
801     h->Sumw2();
802     h->SetXTitle("SPD Clusters");
803     h->SetYTitle("N tracklets");
804     fHistoList->GetList()->Add(h);
805     TH1::AddDirectory(oldStatus);
806   }
807   return h;
808 }
809
810
811 void AliAnalysisTaskTriggerStudy::FillTriggerOverlaps (const char * name, const char * title, Bool_t * vdArray){
812   //Fills a histo with the different trigger statistics in a venn like diagramm. Books it if needed.
813
814   // Get or book histo
815   TString hname = "hTrigStat_";
816   hname+=name;  
817   hname+=fHistoSuffix;
818   TH1 * h = (TH1*) fHistoList->GetList()->FindObject(hname.Data());
819   
820   if(!h) {
821     AliInfo(Form("Booking histo %s",hname.Data()));
822     Bool_t oldStatus = TH1::AddDirectoryStatus();
823     TH1::AddDirectory(kFALSE);
824     Int_t nbins = 0;
825     for(Int_t ientry = 0; ientry < kNVDEntries; ientry++){
826       nbins = nbins | (1<<ientry);
827     }
828     //    cout << "NBINS " << nbins << endl;
829     nbins = nbins+1;
830     h = new TH1I (hname, title, nbins, -0.5, nbins-0.5);
831     fHistoList->GetList()->Add(h);
832     TH1::AddDirectory(oldStatus);
833   
834     // we look at the combinations of n triggers
835     // We set a bit for each trigger to fill the diagram
836     // This is much simpler and faster than any recursive function
837     h->GetXaxis()->SetBinLabel(1,"NONE"); 
838     for(Int_t ibin = 1; ibin < nbins; ibin++){
839       TString binname = "";
840       Bool_t first = kTRUE;
841       for(Int_t ivdentry = 0; ivdentry < kNVDEntries; ivdentry++){
842         if (ibin & (1<<ivdentry)) {
843           if(!first) binname += " & ";
844           binname += kVDNames[ivdentry];
845           first=kFALSE;
846         }
847       }
848       h->GetXaxis()->SetBinLabel(ibin+1,binname.Data());
849     }
850     
851   }
852
853   UInt_t mask = 0;
854   for(Int_t ivdentry = 0; ivdentry < kNVDEntries; ivdentry++){
855     if(vdArray[ivdentry]) {
856       mask  = mask | (1<<ivdentry);
857       //      cout << " 1 "   ;
858     } //else cout << " 0 ";
859   }
860   //  cout << hex << " = " << mask << endl;
861   
862   h->Fill(mask);
863
864 }