Working on the electron cut
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / multPbPb / AliAnalysisTaskMultPbTracks.cxx
1 // AliAnalysisTaskMultPbTracks
2
3 // Author: Michele Floris, CERN
4 // TODO:
5 // - Add chi2/cluster plot for primary, secondaries and fakes
6
7 #include "AliAnalysisTaskMultPbTracks.h"
8 #include "AliESDInputHandler.h"
9 #include "AliAnalysisMultPbTrackHistoManager.h"
10 #include "AliAnalysisManager.h"
11 #include "AliESDtrackCuts.h"
12 #include "AliMCEvent.h"
13 #include "AliStack.h"
14 #include "TH1I.h"
15 #include "TH3D.h"
16 #include "TH2D.h"
17 #include "AliMCParticle.h"
18 #include "AliGenEventHeader.h"
19 #include "AliCentrality.h"
20 #include "AliMultiplicity.h"
21 #include "TFile.h"
22 #include <iostream>
23 #include "AliAnalysisMultPbCentralitySelector.h"
24 #include "AliTriggerAnalysis.h"
25 #include "AliPIDResponse.h"
26
27 using namespace std;
28
29 ClassImp(AliAnalysisTaskMultPbTracks)
30
31 AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks()
32 : AliAnalysisTaskSE("TaskMultPbTracks"),
33   fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fOfflineTrigger(0), fIsMC(0),fIsTPCOnly(0), fTriggerAnalysis(0),fPIDResponse(0),fRejectElectrons(0)
34 {
35   // constructor
36
37   DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
38   DefineOutput(2, AliESDtrackCuts::Class());
39   DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
40
41   fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
42   if(fIsMC) fHistoManager->SetSuffix("MC");
43
44 }
45 AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const char * name)
46   : AliAnalysisTaskSE(name),
47     fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fOfflineTrigger(0),fIsMC(0),fIsTPCOnly(0), fTriggerAnalysis(0),fPIDResponse(0),fRejectElectrons(0)
48 {
49   //
50   // Standard constructur which should be used
51   //
52
53   DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
54   DefineOutput(2, AliESDtrackCuts::Class());
55   DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
56
57   fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
58   if(fIsMC) fHistoManager->SetSuffix("MC");
59
60 }
61
62 AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const AliAnalysisTaskMultPbTracks& obj) : 
63   AliAnalysisTaskSE(obj) ,fESD (0), fHistoManager(0), fCentrSelector(0), fTrackCuts(0),fTrackCutsNoDCA(0),fOfflineTrigger(0),fIsMC(0),fIsTPCOnly(0), fTriggerAnalysis(0),fPIDResponse(0),fRejectElectrons(0)
64 {
65   //copy ctor
66   fESD = obj.fESD ;
67   fHistoManager= obj.fHistoManager;
68   fTrackCuts  = obj.fTrackCuts;
69   fTrackCutsNoDCA  = obj.fTrackCutsNoDCA;
70   fCentrSelector = obj.fCentrSelector;
71   fOfflineTrigger = obj.fOfflineTrigger;
72   fIsMC = obj.fIsMC;
73   fIsTPCOnly = obj.fIsTPCOnly;
74   fTriggerAnalysis = obj.fTriggerAnalysis;
75
76 }
77
78 AliAnalysisTaskMultPbTracks::~AliAnalysisTaskMultPbTracks(){
79   // destructor
80
81   if(!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
82     if(fHistoManager) {
83       delete fHistoManager;
84       fHistoManager = 0;
85       delete fTriggerAnalysis;
86       fTriggerAnalysis=0;
87     }
88   }
89   // Histo list should not be destroyed: fListWrapper is owner!
90
91 }
92 void AliAnalysisTaskMultPbTracks::UserCreateOutputObjects()
93 {
94   // Called once
95
96   // For the DCA distributions, we want to use exactly the same cuts
97   // as for the other distributions, with the exception of the DCA cut
98   fTrackCutsNoDCA = new AliESDtrackCuts(*fTrackCuts); // clone cuts
99   // Reset all DCA cuts; FIXME: is this all?
100   fTrackCutsNoDCA->SetMaxDCAToVertexXY();
101   fTrackCutsNoDCA->SetMaxDCAToVertexZ ();
102   fTrackCutsNoDCA->SetMaxDCAToVertexXYPtDep();
103   fTrackCutsNoDCA->SetMaxDCAToVertexZPtDep();
104
105   fTriggerAnalysis = new AliTriggerAnalysis();
106   fTriggerAnalysis->SetAnalyzeMC(fIsMC);
107
108   //The common PID object can then be retrieved from the input handler. This can naturally be done in the UserCreateOutputObjects:
109   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
110   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
111   fPIDResponse = inputHandler->GetPIDResponse();
112
113   PostData(1,fHistoManager);
114   PostData(2,fTrackCuts);
115   PostData(3,fCentrSelector);
116
117
118 }
119
120
121 void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
122 {
123   // User code
124
125
126   /* PostData(0) is taken care of by AliAnalysisTaskSE */
127   PostData(1,fHistoManager);
128   PostData(2,fTrackCuts);
129   PostData(3,fCentrSelector);
130
131   // Cache (some) histogram pointers to avoid loop in the histo manager
132   static TH3D * hTracks  [AliAnalysisMultPbTrackHistoManager::kNHistos];
133   static TH2D * hDCA     [AliAnalysisMultPbTrackHistoManager::kNHistos];
134   static TH1D * hNTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
135   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
136   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
137   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]    = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
138   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]    = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
139   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]  = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
140   hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
141
142   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
143   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
144   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]    = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
145   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]    = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
146   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]  = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
147   hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
148
149   hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoGen );
150   hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec );
151
152
153   fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Reset();
154
155   fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
156   if (strcmp(fESD->ClassName(),"AliESDEvent")) {
157     AliFatal("Not processing ESDs");
158   }
159   
160   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll);
161
162
163   // Centrality selection
164   Bool_t isCentralitySelected = fCentrSelector->IsCentralityBinSelected(fESD,fTrackCuts);  
165   if(!isCentralitySelected) return;
166
167   // AliESDCentrality *centrality = fESD->GetCentrality();
168   // if(!centrality) {
169   //   AliFatal("Centrality object not available"); 
170   // }
171   // else {
172   //   Int_t centrBin = centrality->GetCentralityClass5(fCentralityEstimator.Data()) ;
173     
174   //   if (centrBin != fCentrBin && fCentrBin != -1) return;
175   // }
176
177   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatCentr);
178
179   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fOfflineTrigger);
180
181   if(!isSelected) return;
182   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatPhysSel);
183
184
185   if (fIsMC) {
186     // Get MC vertex
187     //FIXME: which vertex do I take for MC?
188     TArrayF   vertex;
189     fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
190     Float_t zvGen = vertex[2];    
191
192     if (!fMCEvent) {
193       AliError("No MC info found");
194     } else {
195       
196       //loop on the MC event, only over primaries, which are always
197       //      the first in stack.
198       // WRONG: D&B decays are produced in the transport... Need To Loop over all particles
199       //      Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
200       Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
201       Int_t nPhysicalPrimaries = 0;
202       for (Int_t ipart=0; ipart<nMCTracks; ipart++) { 
203         
204         AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
205         
206         // We don't care about neutrals and non-physical primaries
207         if(mcPart->Charge() == 0) continue;
208
209         // If we are cutting away electrons, also exclude them here:
210         if(fRejectElectrons) 
211           if (TMath::Abs(mcPart->PdgCode())==11) continue;
212
213         //check if current particle is a physical primary and fill map of physical primaries
214         if(!IsPhysicalPrimaryAndTransportBit(ipart)) continue;
215         fHistoManager->FillSpeciesMap(mcPart->PdgCode());
216
217         if(TMath::Abs(mcPart->Zv()-zvGen)>1e-6) {
218           // This cures a bug in Hijing
219           // A little hack here: I put those in the underflow bin of the process type to keep track of them
220           fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(-1);
221           continue;
222         }
223
224         nPhysicalPrimaries++;
225         // Fill species histo and particle species
226         fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID());
227         fHistoManager->FillParticleID(AliAnalysisMultPbTrackHistoManager::kHistoGen, mcPart);
228         
229         //      Float_t zv = vtxESD->GetZ();
230         // Fill generated histo
231         hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zvGen);
232         Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
233         //      cout << "Part " << partCode << endl;
234         if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus  || 
235             partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus ||
236             partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus   || 
237             partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus  ||
238             partCode == AliAnalysisMultPbTrackHistoManager::kPartP       || 
239             partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar  
240             ){
241           //      cout << "Found " << partCode << endl;
242           
243           fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen, partCode);
244         }
245       }
246       hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries);
247       fHistoManager->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(zvGen);
248   
249     
250     }
251   }
252   
253
254   const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
255   if(!vtxESD) return;
256   // TODO: leave the cuts here or find a better place?)
257   // Quality cut on vertexer Z, as suggested by Francesco Prino
258   if(vtxESD->IsFromVertexerZ()) {
259     if (vtxESD->GetNContributors() <= 0) return;
260     if (vtxESD->GetDispersion() > 0.04) return;
261     if (vtxESD->GetZRes() > 0.25) return;
262   }
263   // "Beam gas" vertex cut
264   const AliESDVertex * vtxESDTPC= fESD->GetPrimaryVertexTPC(); 
265   if(vtxESDTPC->GetNContributors()<1) return;
266   if (vtxESDTPC->GetNContributors()<(-10.+0.25*fESD->GetMultiplicity()->GetNumberOfITSClusters(0)))     return;  
267
268   // Fill  statistics
269   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtx);
270   
271   // ZDC cut, only ZNs
272   Bool_t zdcA   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kASide, kTRUE, kFALSE) ; 
273   Bool_t zdcC   = fTriggerAnalysis->ZDCTDCTrigger(fESD, AliTriggerAnalysis::kCSide, kTRUE, kFALSE) ;                    
274
275   if (!(zdcA && zdcC) && (!fIsMC)) return;
276   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatZDCCut);
277
278
279   if(TMath::Abs(vtxESD->GetZ()) > 10) return;
280   fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtxRangeCut);
281
282   // FIXME
283   Float_t ntracksOk = fTrackCuts->CountAcceptedTracks(fESD);
284   const AliMultiplicity * mult = fESD->GetMultiplicity();
285   if (ntracksOk < 1000 &&  mult->GetNumberOfITSClusters(1) > 4500) {
286     Float_t dummy;
287     Printf("IEV: %d, Orbit: %x, Period: %d, BC: %d\n",fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber());
288     printf("%s ----> Processing event # %lld\n", CurrentFileName(), Entry());
289     cout << "File " << AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName() << endl;   
290     cout << "Nt " << ntracksOk << ", "<< fESD->GetNumberOfTracks() << " V0 " <<  fCentrSelector->GetCorrV0(fESD,dummy)  << " SPD1 " << mult->GetNumberOfITSClusters(1) << endl;       
291   }
292
293
294   
295   // Fill Vertex
296
297   fHistoManager->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(vtxESD->GetZ());
298
299   // loop on tracks
300   Int_t ntracks = fESD->GetNumberOfTracks();
301   Int_t acceptedTracks = 0;
302
303   for (Int_t itrack = 0; itrack<ntracks; itrack++) {    
304     AliESDtrack * esdTrack = fIsTPCOnly ? fTrackCuts->GetTPCOnlyTrack(fESD,itrack) : fESD->GetTrack(itrack);// FIXMEL TrackCuts or TrackCutsNoDCA for the TPC?
305     if (!esdTrack) continue;
306           
307     // Fill DCA distibutions, without the DCA cut, Fill the other stuff, with all cuts!
308     Bool_t accepted = fTrackCuts->AcceptTrack(esdTrack);
309     Bool_t acceptedNoDCA = fTrackCutsNoDCA->AcceptTrack(esdTrack);
310
311     //fHistoManager->GetHistoElectronCutQA()->Fill(esdTrack->Pt(), esdTrack->GetTPCsignal()); FIXME: temporary hack
312
313     if(fRejectElectrons && acceptedNoDCA) { // don't check this for tracks to be anyways rejected
314       Bool_t isEl = IsElectron(esdTrack);
315       if(isEl) {
316         //      cout << " --> FOUND!" << endl;  
317         fHistoManager->GetHistoElectronCutQA()->Fill(esdTrack->Pt(), esdTrack->GetTPCsignal());
318       }
319       accepted = accepted && !isEl;
320       acceptedNoDCA = acceptedNoDCA && !isEl;
321     }
322
323     if(accepted) acceptedTracks++;
324
325     // Compute weighted offset
326     // TODO: other choiches for the DCA?
327     Double_t weightedDCA = 10;
328     
329
330 #if defined WEIGHTED_DCA
331     Double_t b = fESD->GetMagneticField();
332     Double_t dca[2];
333     Double_t cov[3];
334     if (esdTrack->PropagateToDCA(vtxESD, b,3., dca, cov)) {
335       Double_t det = cov[0]*cov[2]-cov[1]*cov[1]; 
336       if (det<=0) {
337         AliError("DCA Covariance matrix is not positive definite!");
338       }
339       else {
340         weightedDCA = (dca[0]*dca[0]*cov[2]+dca[1]*dca[1]*cov[0]-2*dca[0]*dca[1]*cov[1])/det; 
341         weightedDCA = weightedDCA>0 ? TMath::Sqrt(weightedDCA) : 10;
342       }
343       //      printf("dR:%e dZ%e  sigR2:%e sigRZ:%e sigZ2:%e\n",dca[0],dca[1],cov[0],cov[1],cov[2]);
344     }
345     
346 #elif defined TRANSVERSE_DCA
347     Float_t xz[2]; 
348     esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz); 
349     weightedDCA = xz[0];
350 #endif
351     
352
353     // Alternative: p*DCA
354     // Float_t xz[2]; 
355     // esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz); 
356     // Float_t dca = xz[0]*esdTrack->P();
357
358   // This cures a bug in Hijing (displaced primaries)
359   if (fIsMC) {
360     Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!    
361     if (IsPhysicalPrimaryAndTransportBit(label)) {
362       AliMCParticle *mcPart  =  (AliMCParticle*)fMCEvent->GetTrack(label);
363       if(!mcPart) AliFatal("No particle");
364       TArrayF   vertex;
365       fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
366       Float_t zvGen = vertex[2];    
367       if(TMath::Abs(mcPart->Zv()-zvGen)>1e-6) continue;    
368     }
369   }
370   // for each track
371   if(accepted) {
372     hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
373     if (TMath::Abs(esdTrack->Eta())<0.5 && TMath::Abs(vtxESD->GetZ())<7)
374       fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(esdTrack->Pt());
375   }
376   if(acceptedNoDCA)
377     hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA,esdTrack->Pt());
378
379     // Fill separately primaries and secondaries
380     // FIXME: fakes? Is this the correct way to account for them?
381     // Get label and corresponding mcPart;
382     if (fIsMC) {
383       Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
384       AliMCParticle *mcPart  = label < 0 ? 0 : (AliMCParticle*)fMCEvent->GetTrack(label);
385       if (!mcPart)  {
386         if(accepted)
387           hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
388         if(acceptedNoDCA)
389           hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA,esdTrack->Pt());
390       }
391       else {
392         if(IsPhysicalPrimaryAndTransportBit(label)) {
393           // Fill species histo
394           fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->Fill(mcPart->Particle()->GetUniqueID());
395           if(accepted) {
396             hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
397             Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
398             if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus  || 
399                 partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus ||
400                 partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus   || 
401                 partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus  ||
402                 partCode == AliAnalysisMultPbTrackHistoManager::kPartP       || 
403                 partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar    ||
404                 partCode == AliAnalysisMultPbTrackHistoManager::kPartLPlus   || 
405                 partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus 
406                 )
407               fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, partCode);
408           }
409           if(acceptedNoDCA)
410             hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(weightedDCA,esdTrack->Pt());
411         } 
412         else {
413           Int_t mfl=0;
414           Int_t indexMoth=mcPart->Particle()->GetFirstMother();
415           if(indexMoth>=0){
416             TParticle* moth = fMCEvent->Stack()->Particle(indexMoth);
417             Float_t codemoth = TMath::Abs(moth->GetPdgCode());
418             mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
419           }
420           if(mfl==3){ // strangeness
421             fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak)->Fill(mcPart->Particle()->GetUniqueID());
422             if(accepted)
423               hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
424             if(acceptedNoDCA)
425               hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA,esdTrack->Pt());       
426           }else{ // material
427             fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat)->Fill(mcPart->Particle()->GetUniqueID());
428             if(accepted)
429               hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
430             if(acceptedNoDCA)
431               hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(weightedDCA,esdTrack->Pt());        
432
433           }
434
435
436         }
437       }
438     }
439
440
441   }
442   //  cout << acceptedTracks << endl;
443   
444   hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]  ->Fill(acceptedTracks);
445
446   // FIXME: this used to be filled with rescaled V0. I'm using raw V0 to test some newer production here.
447   AliESDVZERO* esdV0 = fESD->GetVZEROData();
448   Float_t multV0A=esdV0->GetMTotV0A();
449   Float_t multV0C=esdV0->GetMTotV0C();
450   Float_t multV0 = multV0A+multV0C;
451   fHistoManager->GetHistoV0vsNtracks(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(acceptedTracks, multV0);
452   // FIXME: the old version:
453   // Float_t v0;
454   // fHistoManager->GetHistoV0vsNtracks(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(acceptedTracks, fCentrSelector->GetCorrV0(fESD,v0));
455
456   // Fill <pt> analysis histos
457   fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Scale(1,"width");// correct bin width
458   if (fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetEntries()>0)
459     fHistoManager->GetHistoMeanPt(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetMean());
460   if (fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetMean() >  
461       fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecHighestMeanPt)->GetMean()) {
462     // Found a new highest pt: first reset than add it to the container for the highest
463     fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecHighestMeanPt)->Reset();
464     fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecHighestMeanPt)
465       ->Add(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec));
466   }
467   if ((fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetMean() <  
468        fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)->GetMean() &&
469        fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->GetEntries()>0) ||
470       !(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)->GetEntries()>0)) {
471     // Found a new lowest pt: first reset than add it to the container for the lowest
472     fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)->Reset();
473     fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRecLowestMeanPt)
474       ->Add(fHistoManager->GetHistoPtEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec));
475   }
476 }
477
478 void   AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
479   // terminate
480
481 }
482
483
484 Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) {
485   // Check if it is a primary and if it was transported
486   Bool_t physprim=fMCEvent->IsPhysicalPrimary(ipart);
487   if (!physprim) return kFALSE;
488   Bool_t transported = ((AliMCParticle*)fMCEvent->GetTrack(ipart))->Particle()->TestBit(kTransportBit);
489   if(!transported) return kFALSE;
490
491   return kTRUE;
492
493 }
494
495 Bool_t AliAnalysisTaskMultPbTracks::IsElectron(AliESDtrack * esdTrack) {
496   Bool_t isElectron = kFALSE;
497     // ((fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kElectron) < 2 ) &&
498     //  (fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kPion    ) > 1 ) && 
499     //  (fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kProton  ) > 1 ) && 
500     //  (fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kKaon    ) > 1 ) 
501     //  );//FIXME SKIP ELECTRONS below p = 1.2 gev, make configurable. Keep particles if they are in the crossing
502
503   // cout << " - TPC: " << fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kElectron) << ", TOF: " << fPIDResponse->NumberOfSigmasTOF(esdTrack,AliPID::kElectron);
504   // printf("[%f,%f,%d]\n", esdTrack->Pt(), esdTrack->Eta(), esdTrack->GetTPCclusters(0));  
505   // if (fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kElectron) <0){
506
507     
508   // } 
509   isElectron = (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(esdTrack,AliPID::kElectron)) < 2) && (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(esdTrack,AliPID::kElectron) < 3));
510
511   return isElectron;
512 }
513
514 // void AliAnalysisTaskEvil::PrintProcInfo()
515 // {
516 //   ProcInfo_t info;
517 //   gSystem->GetProcInfo(&info);
518 //   AliInfo(Form("fMemResident=%10ld kB  fMemVirtual=%10ld kB",info.fMemResident,info.fMemVirtual));
519 // }