rule checker simple fixes and cosmetics
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / AliAnalysisTaskCounter.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //_________________________________________________________________________
17 // Count events with different selections
18 //
19 // It produces a histogram with the number of events with 9 bins:
20 // 0: all events (that passed the physics selection if it was on)
21 // 1: same but cross check that event pointer did exist
22 // 2: passes vertex cut
23 // 3: passes track number cut, tracks for eta < 0.8
24 // 4: 3 && 2
25 // 5: pass VAND
26 // 6: 5 && 2
27 // 7: 5 && 3
28 // 8: 5 && 3 && 2
29 // 9: not pileup from SPD
30 // 10: Good vertex
31 // 11: 10 && 5
32 // 12: 10 && 3
33 // 13: 10 && 2
34 // 14: 10 && 2 && 3 && 5
35 // 15: 10 && 9
36 // 16: 9  && 5
37 //
38 // Author: Gustavo Conesa Balbastre (LPSC)
39 //         
40 //_________________________________________________________________________
41
42 #include "TH2F.h"
43 #include "AliAODHeader.h"
44 #include "AliTriggerAnalysis.h"
45 #include "AliESDEvent.h"
46 #include "AliAODEvent.h"
47 #include "AliESDtrackCuts.h"
48 #include "AliAnalysisManager.h"
49 #include "AliInputEventHandler.h"
50
51 #include "AliAnalysisTaskCounter.h"
52 ClassImp(AliAnalysisTaskCounter)
53
54 //______________________________________________________________
55 AliAnalysisTaskCounter::AliAnalysisTaskCounter(const char *name) 
56 : AliAnalysisTaskSE(name), 
57   fAcceptFastCluster(kTRUE),
58   fZVertexCut(10.), 
59   fTrackMultEtaCut(0.8),
60   fCaloFilterPatch(kFALSE),
61   fOutputContainer(0x0), 
62   fESDtrackCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
63   fTriggerAnalysis (new AliTriggerAnalysis),
64   fhNEvents(0),
65   fhXVertex(0),    fhYVertex(0),    fhZVertex(0),
66   fhXGoodVertex(0),fhYGoodVertex(0),fhZGoodVertex(0)
67 {
68   //ctor
69   DefineOutput(1, TList::Class());
70 }
71
72 //______________________________________________
73 AliAnalysisTaskCounter::AliAnalysisTaskCounter() 
74   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskCounter"),
75     fAcceptFastCluster(kTRUE),
76     fZVertexCut(10.),
77     fTrackMultEtaCut(0.8),
78     fCaloFilterPatch(kFALSE),
79     fOutputContainer(0x0), 
80     fESDtrackCuts(AliESDtrackCuts::GetStandardITSTPCTrackCuts2010()),
81     fTriggerAnalysis (new AliTriggerAnalysis),
82     fhNEvents(0),
83     fhXVertex(0),    fhYVertex(0),    fhZVertex(0),
84     fhXGoodVertex(0),fhYGoodVertex(0),fhZGoodVertex(0)
85 {
86   // ctor
87   DefineOutput(1, TList::Class());
88 }
89
90 //__________________________________________________
91 AliAnalysisTaskCounter::~AliAnalysisTaskCounter()
92 {
93   //Destructor
94   if(fOutputContainer){
95     fOutputContainer->Delete() ; 
96     delete fOutputContainer ;
97   }
98   
99   if(fESDtrackCuts)    delete fESDtrackCuts;
100   if(fTriggerAnalysis) delete fTriggerAnalysis;
101   
102 }
103
104
105 //____________________________________________________
106 void AliAnalysisTaskCounter::UserCreateOutputObjects()
107 {
108   // Init histograms
109   
110   fOutputContainer = new TList();
111   
112   fhZVertex     = new TH1F("hZVertex", " Z vertex distribution"   , 200 , -50 , 50  ) ;
113   fhZVertex->SetXTitle("v_{z} (cm)");
114   fOutputContainer->Add(fhZVertex);
115
116   fhZGoodVertex     = new TH1F("hZGoodVertex", " Good Z vertex distribution"   , 200 , -50 , 50  ) ;
117   fhZGoodVertex->SetXTitle("v_{z} (cm)");
118   fOutputContainer->Add(fhZGoodVertex);
119   
120   fhXVertex     = new TH1F("hXVertex", " X vertex distribution"   , 200 , -2 , 2  ) ;
121   fhXVertex->SetXTitle("v_{x} (cm)");
122   fOutputContainer->Add(fhXVertex);
123   
124   fhXGoodVertex     = new TH1F("hXGoodVertex", " Good X vertex distribution"   , 200 , -2 , 2  ) ;
125   fhXGoodVertex->SetXTitle("v_{x} (cm)");
126   fOutputContainer->Add(fhXGoodVertex);
127   
128   fhYVertex     = new TH1F("hYVertex", " Y vertex distribution"   , 200 , -2 , 2  ) ;
129   fhYVertex->SetXTitle("v_{y} (cm)");
130   fOutputContainer->Add(fhYVertex);
131   
132   fhYGoodVertex     = new TH1F("hYGoodVertex", " Good Y vertex distribution"   , 200 , -2 , 2  ) ;
133   fhYGoodVertex->SetXTitle("v_{y} (cm)");
134   fOutputContainer->Add(fhYGoodVertex);
135   
136   
137   fhNEvents = new TH1I("hNEvents", "Number of analyzed events", 21, 0, 21) ;
138   fhNEvents->SetXTitle("Selection");
139   fhNEvents->SetYTitle("# events");
140   fhNEvents->GetXaxis()->SetBinLabel(1 ,"1  = PS");
141   fhNEvents->GetXaxis()->SetBinLabel(2 ,"2  = 1  & ESD");
142   fhNEvents->GetXaxis()->SetBinLabel(3 ,"3  = 2  & |Z|<10");
143   fhNEvents->GetXaxis()->SetBinLabel(4 ,"4  = 2  & |eta|<0.8");
144   fhNEvents->GetXaxis()->SetBinLabel(5 ,"5  = 3  & 4");
145   fhNEvents->GetXaxis()->SetBinLabel(6 ,"6  = 2  & V0AND");
146   fhNEvents->GetXaxis()->SetBinLabel(7 ,"7  = 3  & 6");
147   fhNEvents->GetXaxis()->SetBinLabel(8 ,"8  = 4  & 6");
148   fhNEvents->GetXaxis()->SetBinLabel(9 ,"9  = 5  & 6");
149   fhNEvents->GetXaxis()->SetBinLabel(10,"10 = 2  & not pileup");
150   fhNEvents->GetXaxis()->SetBinLabel(11,"11 = 2  & good vertex");
151   fhNEvents->GetXaxis()->SetBinLabel(12,"12 = 3  & 11");
152   fhNEvents->GetXaxis()->SetBinLabel(13,"13 = 4  & 11");
153   fhNEvents->GetXaxis()->SetBinLabel(14,"14 = 6  & 11");
154   fhNEvents->GetXaxis()->SetBinLabel(15,"15 = 9  & 11");
155   fhNEvents->GetXaxis()->SetBinLabel(16,"16 = 10 & 11");
156   fhNEvents->GetXaxis()->SetBinLabel(17,"17 = 6  & 10");
157   fhNEvents->GetXaxis()->SetBinLabel(17,"17 = 1  & |Z|<50");  
158   fhNEvents->GetXaxis()->SetBinLabel(18,"18 = Reject EMCAL 1");
159   fhNEvents->GetXaxis()->SetBinLabel(19,"19 = 18 & 2");
160   fhNEvents->GetXaxis()->SetBinLabel(20,"20 = Reject EMCAL 2");
161   fhNEvents->GetXaxis()->SetBinLabel(21,"20 = 20 & 2");
162
163   fOutputContainer->Add(fhNEvents);
164
165   fOutputContainer->SetOwner(kTRUE);
166   
167   PostData(1,fOutputContainer);
168   
169 }
170
171 //_______________________________________________
172 void AliAnalysisTaskCounter::UserExec(Option_t *) 
173 {
174   // Main loop
175   // Called for each event
176   
177   //printf("___ Event __ %d __\n",(Int_t)Entry());
178   
179   fhNEvents->Fill(0.5);  
180   
181   AliVEvent * event = InputEvent();
182   if (!event) 
183   {
184     printf("AliAnalysisTaskCounter::UserExec() - ERROR: event not available \n");
185     return;
186   }
187   
188   AliESDEvent * esdevent = dynamic_cast<AliESDEvent*> (event);
189   AliAODEvent * aodevent = dynamic_cast<AliAODEvent*> (event);
190   
191   TString triggerclasses = "";
192   if(esdevent) triggerclasses = esdevent->GetFiredTriggerClasses();
193   if(aodevent) triggerclasses = aodevent->GetFiredTriggerClasses();
194
195   if (triggerclasses.Contains("FAST") && !triggerclasses.Contains("ALL") && !fAcceptFastCluster) 
196   {
197     //printf("Do not count events from fast cluster, trigger name %s\n",triggerclasses.Data());
198     return;
199   }
200
201   fhNEvents->Fill(1.5);  
202     
203   //Initialize bools
204   Bool_t bSelectVZ    = kFALSE;
205   Bool_t bV0AND       = kFALSE; 
206   Bool_t bPileup      = kFALSE;
207   Bool_t bGoodV       = kFALSE;
208   Bool_t bSelectTrack = kFALSE;   
209   Int_t  trackMult    = 0;
210   
211   //---------------------------------
212   //Get the primary vertex, cut on Z
213   //---------------------------------
214   Double_t v[3];
215   event->GetPrimaryVertex()->GetXYZ(v) ;
216   fhXVertex->Fill(v[0]);
217   fhYVertex->Fill(v[1]);
218   fhZVertex->Fill(v[2]);
219   
220   if(TMath::Abs(v[2]) < fZVertexCut) 
221   {
222     bSelectVZ=kTRUE;
223     fhNEvents->Fill(2.5);  
224   }
225   //else printf("Vertex out %f \n",v[2]);
226   
227
228   //--------------------------------------------------
229   //Tweak for calorimeter only productions
230   //--------------------------------------------------
231   if(fCaloFilterPatch && !esdevent)
232   { 
233     if(event->GetNumberOfCaloClusters() > 0) 
234     {
235       AliVCluster * calo = event->GetCaloCluster(0);
236       if(calo->GetNLabels() == 4){
237         Int_t * selection = calo->GetLabels();
238         bPileup   = selection[0];
239         bGoodV    = selection[1]; 
240         bV0AND    = selection[2]; 
241         trackMult = selection[3];
242         //if(selection[0] || selection[1] || selection[2])
243         //printf(" pu %d, gv %d, v0 %d, track mult %d\n ", selection[0], selection[1], selection[2], selection[3]);
244         if(trackMult > 0 )  
245           bSelectTrack = kFALSE;
246       } 
247       else 
248       {
249         //First filtered AODs, track multiplicity stored there.  
250         trackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
251       }
252     }
253     else
254     {   //at least one cluster
255         //printf("AliAnalysisTaskCounter::UserExec() - No clusters in event\n");
256         //Remove events with  vertex (0,0,0), bad vertex reconstruction
257         if(TMath::Abs(v[0]) < 1.e-6 && TMath::Abs(v[1]) < 1.e-6 && TMath::Abs(v[2]) < 1.e-6) bGoodV = kFALSE;
258       
259         //First filtered AODs, track multiplicity stored there.  
260         trackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
261     }
262   }
263   else 
264   {
265     //--------------------------------------------------
266     //Count tracks, cut on number of tracks in eta < 0.8
267     //--------------------------------------------------
268     Int_t nTracks   = event->GetNumberOfTracks() ;
269     for (Int_t itrack =  0; itrack <  nTracks; itrack++) 
270     {////////////// track loop
271       AliVTrack * track = (AliVTrack*)event->GetTrack(itrack) ; // retrieve track from esd
272       
273       //Only for ESDs
274       if(esdevent && !fESDtrackCuts->AcceptTrack((AliESDtrack*)track)) continue;
275       
276       //Do not count tracks out of acceptance cut
277       if(TMath::Abs(track->Eta())< fTrackMultEtaCut) trackMult++;
278     }
279   }
280   
281   //printf("AliAnalysisTaskCounter::UserExec() - Track Mult %d \n",trackMult);
282   
283   //--------------------------------------------------
284   // At least one track
285   //--------------------------------------------------
286   if (trackMult > 0) {
287     bSelectTrack = kTRUE; 
288                   fhNEvents->Fill(3.5);
289     if(bSelectVZ) fhNEvents->Fill(4.5);
290   }
291   
292   //---------------------------------
293   // V0AND
294   //---------------------------------
295   if(esdevent && !fCaloFilterPatch) bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esdevent, AliTriggerAnalysis::kV0AND);
296   //else if(aodevent  && !fCaloFilterPatch) bV0AND = //FIXME FOR AODs
297   
298   if(bV0AND)
299   {
300                                     fhNEvents->Fill(5.5);
301     if (bSelectVZ)                  fhNEvents->Fill(6.5);
302     if (bSelectTrack)               fhNEvents->Fill(7.5);
303     if (bSelectVZ &&  bSelectTrack) fhNEvents->Fill(8.5);
304   }
305   
306   //---------------------------------
307   // Pileup
308   //---------------------------------
309   if(!fCaloFilterPatch)
310     bPileup = event->IsPileupFromSPD(3, 0.8, 3., 2., 5.); //Default values, if not it does not compile
311   //bPileup = event->IsPileupFromSPD(); 
312   
313   if (!bPileup)
314   {
315                 fhNEvents->Fill(9.5);
316     if(bV0AND)  fhNEvents->Fill(16.5);
317   }
318   
319   //---------------------------------
320   // Good vertex
321   //---------------------------------
322   if(esdevent && !fCaloFilterPatch) bGoodV = CheckForPrimaryVertex();
323   if(bGoodV) 
324   {
325     fhXGoodVertex->Fill(v[0]);
326     fhYGoodVertex->Fill(v[1]);
327     fhZGoodVertex->Fill(v[2]);
328     
329                      fhNEvents->Fill(10.5);
330     if(bSelectVZ)    fhNEvents->Fill(11.5);
331     if(bSelectTrack) fhNEvents->Fill(12.5);
332     if(bV0AND)       fhNEvents->Fill(13.5);
333     if(bSelectVZ && bSelectTrack && bV0AND)    
334                      fhNEvents->Fill(14.5); 
335     if(!bPileup)     fhNEvents->Fill(15.5); 
336   }
337
338   //printf("AliAnalysisTaskCounter::UserExec() : z vertex %d, good vertex %d, v0and %d, pile up %d, track mult %d\n ", bSelectVZ, bGoodV, bV0AND, bPileup, trackMult);
339   
340   // Events that could be rejected in EMCAL
341   // LHC11a, SM4 and some SM3 events cut with this
342   Bool_t bEMCALRejected = kFALSE;
343   for (Int_t i = 0; i < InputEvent()->GetNumberOfCaloClusters(); i++)
344   {
345     AliVCluster *clus = InputEvent()->GetCaloCluster(i);
346     if(clus->IsEMCAL()){    
347       if ((clus->E() > 500 && clus->GetNCells() > 200 ) || clus->GetNCells() > 200)  
348       {
349         
350         //printf("Counter: Reject event with cluster: E %f, ncells %d\n",clus->E(),clus->GetNCells());
351         
352                          fhNEvents->Fill(17.5); 
353         if(bSelectVZ)    fhNEvents->Fill(18.5);
354         bEMCALRejected = kTRUE;
355         break;
356       }
357     }
358   }
359   
360   //LHC11a, 3 last runs, cut with this
361   if(!bEMCALRejected)
362   {
363     // Count number of cells in SM3 with energy larger than 0.1, cut on this number
364     Int_t ncellsSM3 = 0;
365     Int_t ncellsSM4 = 0;
366     for(Int_t icell = 0; icell < event->GetEMCALCells()->GetNumberOfCells(); icell++)
367     {
368       if(event->GetEMCALCells()->GetAmplitude(icell) > 0.1 && event->GetEMCALCells()->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
369       if(event->GetEMCALCells()->GetAmplitude(icell) > 0.1 && event->GetEMCALCells()->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;
370     }
371     
372     Int_t ncellcut = 21;
373     if(triggerclasses.Contains("EMC")) ncellcut = 35;
374     
375     if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 )
376     {
377       //printf("Counter: reject event with ncells in SM3: ncells %d\n",ncells);
378
379                        fhNEvents->Fill(19.5); 
380       if(bSelectVZ)    fhNEvents->Fill(20.5);
381     }
382     
383   }
384   
385   PostData(1,fOutputContainer);
386
387 }
388
389 //____________________________________________________
390 Bool_t AliAnalysisTaskCounter::CheckForPrimaryVertex()
391 {
392   //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
393   //It only works for ESDs
394   
395   AliESDEvent * event = dynamic_cast<AliESDEvent*> (InputEvent());
396   if(!event) return 0;
397   
398   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0) 
399   {
400     return 1;
401   }
402   
403   if(event->GetPrimaryVertexTracks()->GetNContributors() < 1) 
404   {
405     // SPD vertex
406     if(event->GetPrimaryVertexSPD()->GetNContributors() > 0) 
407     {
408       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
409       return 1;
410       
411     }
412     if(event->GetPrimaryVertexSPD()->GetNContributors() < 1) 
413     {
414       //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
415       return 0;
416     }
417   }
418   
419   return 0;
420   //return fInputEvent->GetPrimaryVertex()->GetNContributors()>0;
421 }
422
423
424
425 //_____________________________________________________
426 void AliAnalysisTaskCounter::FinishTaskOutput()
427 {
428   // Put in the output some event summary histograms
429   
430   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
431   AliInputEventHandler *inputH = dynamic_cast<AliInputEventHandler*>(am->GetInputEventHandler());
432   if (!inputH) return; 
433   TH2F *histStat = dynamic_cast<TH2F*>(inputH->GetStatistics()); 
434   TH2F *histBin0 = dynamic_cast<TH2F*>(inputH->GetStatistics("BIN0"));
435   
436   if(histStat)
437     fOutputContainer->Add(histStat);
438   else if(DebugLevel() > 1)
439     printf("AliAnalysisTaskCounter::FinishTaskOutput() - Stat histogram not available check, \n if ESDs, that AliPhysicsSelection was on, \n if AODs, if EventStat_temp.root exists \n");
440
441   if(histBin0)
442     fOutputContainer->Add(histBin0); 
443   
444 }