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