update to AOD analysis
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskSpectraAOD.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 \r
16 //-----------------------------------------------------------------\r
17 //         AliAnalysisTaskSpectraAOD class\r
18 //-----------------------------------------------------------------\r
19 \r
20 #include "TChain.h"\r
21 #include "TTree.h"\r
22 #include "TLegend.h"\r
23 #include "TH1F.h"\r
24 #include "TH2F.h"\r
25 #include "TCanvas.h"\r
26 #include "AliAnalysisTask.h"\r
27 #include "AliAnalysisManager.h"\r
28 #include "AliAODTrack.h"\r
29 #include "AliAODMCParticle.h"\r
30 #include "AliAODEvent.h"\r
31 #include "AliAODInputHandler.h"\r
32 #include "AliAnalysisTaskSpectraAOD.h"\r
33 #include "AliAnalysisTaskESDfilter.h"\r
34 #include "AliAnalysisDataContainer.h"\r
35 #include "AliSpectraAODHistoManager.h"\r
36 #include "AliSpectraAODTrackCuts.h"\r
37 #include "AliSpectraAODEventCuts.h"\r
38 #include "AliCentrality.h"\r
39 #include "TProof.h"\r
40 #include "AliPID.h"\r
41 #include "AliVEvent.h"\r
42 #include "AliPIDResponse.h"\r
43 #include "AliStack.h"\r
44 #include <iostream>\r
45 \r
46 \r
47 \r
48 \r
49 using namespace AliSpectraNameSpace;\r
50 using namespace std;\r
51 \r
52 ClassImp(AliAnalysisTaskSpectraAOD) // EX1 This stuff tells root to implement the streamer, inspector methods etc (we discussed about it today)\r
53 \r
54 //________________________________________________________________________\r
55 AliAnalysisTaskSpectraAOD::AliAnalysisTaskSpectraAOD(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0), fIsMC(0), fPIDResponse(0), fNSigmaPID(0), fYCut(0)\r
56 {\r
57   // Default constructor\r
58   fNSigmaPID = 3.;\r
59   fYCut = .5;\r
60   DefineInput(0, TChain::Class());\r
61   DefineOutput(1, AliSpectraAODHistoManager::Class());\r
62   DefineOutput(2, AliSpectraAODEventCuts::Class());\r
63   DefineOutput(3, AliSpectraAODTrackCuts::Class());\r
64 \r
65 }\r
66 //________________________________________________________________________\r
67 Bool_t AliAnalysisTaskSpectraAOD::CheckYCut(AODParticleSpecies_t species, AliAODTrack* track) const\r
68 {\r
69   // check if the rapidity is within the set range\r
70   // note: masses are hardcoded for now. we could look them up in the pdg database, but that would mean accecing it 100k+ times per run ...\r
71   Double_t y;\r
72   if (species == kProton) { y = track->Y(9.38271999999999995e-01); }\r
73   if ( species == kKaon ) { y = track->Y(4.93676999999999977e-01); }\r
74   if ( species == kPion)  { y = track->Y(1.39570000000000000e-01); }\r
75   if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;\r
76   return kTRUE;\r
77 }\r
78 //____________________________________________________________________________\r
79 Bool_t AliAnalysisTaskSpectraAOD::CheckYCut(AliAODMCParticle* particle) const\r
80 {\r
81   // check if the rapidity is within the set range\r
82   Double_t y = particle->Y();\r
83   if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;\r
84   return kTRUE;\r
85 }\r
86 //________________________________________________________________________\r
87 void AliAnalysisTaskSpectraAOD::UserCreateOutputObjects()\r
88 {\r
89   // create output objects\r
90   fHistMan = new AliSpectraAODHistoManager("SpectraHistos");\r
91 \r
92   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
93   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
94 \r
95   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();\r
96   AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());\r
97   fPIDResponse = inputHandler->GetPIDResponse();\r
98 \r
99   PostData(1, fHistMan);\r
100   PostData(2, fEventCuts);\r
101   PostData(3, fTrackCuts);\r
102 \r
103 }\r
104 //________________________________________________________________________\r
105 void AliAnalysisTaskSpectraAOD::UserExec(Option_t *)\r
106 {\r
107   // main event loop\r
108   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);\r
109   if (strcmp(fAOD->ClassName(), "AliAODEvent"))\r
110     {\r
111       AliFatal("Not processing AODs");\r
112     }\r
113   \r
114   if(!fPIDResponse) {\r
115     AliError("Cannot get pid response");\r
116     return;\r
117   }\r
118   \r
119   //Selection on QVector, before ANY other selection on the event\r
120   //Spectra MUST be normalized wrt events AFTER the selection on Qvector\r
121   // Can we include this in fEventCuts\r
122   Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;\r
123   Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;\r
124   Int_t multPos = 0;\r
125   Int_t multNeg = 0;\r
126   for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {\r
127     AliAODTrack* aodTrack = fAOD->GetTrack(iT);\r
128     if (!fTrackCuts->IsSelected(aodTrack)) continue;\r
129     if (aodTrack->Eta() >= 0){\r
130       multPos++;\r
131       Qx2EtaPos += TMath::Cos(2*aodTrack->Phi()); \r
132       Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());\r
133     } else {\r
134       multNeg++;\r
135       Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi()); \r
136       Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());\r
137     }\r
138   } \r
139   Double_t qPos = TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);\r
140   Double_t qNeg = TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);\r
141   if((qPos>fTrackCuts->GetQvecMin() && qPos<fTrackCuts->GetQvecMax()) || (qNeg>fTrackCuts->GetQvecMin() && qNeg<fTrackCuts->GetQvecMax())){\r
142     \r
143     //check on centrality distribution\r
144     fHistMan->GetPtHistogram("CentCheck")->Fill(fAOD->GetCentrality()->GetCentralityPercentile("V0M"),fAOD->GetHeader()->GetCentralityP()->GetCentralityPercentileUnchecked("V0M"));\r
145     \r
146     if (!fEventCuts->IsSelected(fAOD)) return;\r
147     \r
148     //fill q distributions vs centrality, after all event selection\r
149     fHistMan->GetqVecHistogram(kHistqVecPos)->Fill(qPos,fAOD->GetCentrality()->GetCentralityPercentile("V0M"));  // qVector distribution\r
150     fHistMan->GetqVecHistogram(kHistqVecNeg)->Fill(qNeg,fAOD->GetCentrality()->GetCentralityPercentile("V0M"));  // qVector distribution\r
151     \r
152     //AliCentrality fAliCentral*;\r
153     //   if ((fAOD->GetCentrality())->GetCentralityPercentile("V0M") > 5.) return;\r
154     \r
155     // First do MC to fill up the MC particle array, such that we can use it later\r
156     TClonesArray *arrayMC = 0;\r
157     if (fIsMC)\r
158       {\r
159         arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
160         if (!arrayMC) {\r
161           AliFatal("Error: MC particles branch not found!\n");\r
162         }\r
163         Int_t nMC = arrayMC->GetEntries();\r
164         for (Int_t iMC = 0; iMC < nMC; iMC++)\r
165           {\r
166             AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
167             fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
168             \r
169             if(TMath::Abs(partMC->Eta()) > fTrackCuts->GetEta()) continue;\r
170             // check for true PID + and fill P_t histos \r
171             //if (partMC->IsPhysicalPrimary() && CheckYCut(partMC) ) {// only primary vertices and y cut satisfied\r
172             if (CheckYCut(partMC) ){// only primary vertices and y cut satisfied\r
173               if ( partMC->PdgCode() == 2212) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryProtonPlus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
174               if ( partMC->PdgCode() == -2212) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryProtonMinus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
175               if ( partMC->PdgCode() == 321)  { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryKaonPlus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
176               if ( partMC->PdgCode() == -321) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryKaonMinus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
177               if ( partMC->PdgCode() == 211)  { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryPionPlus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
178               if ( partMC->PdgCode() == -211) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryPionMinus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
179             }\r
180           }\r
181       }\r
182     \r
183     \r
184     \r
185     //main loop on tracks\r
186     for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)\r
187       {\r
188         \r
189         AliAODTrack* track = fAOD->GetTrack(iTracks);\r
190         if (!fTrackCuts->IsSelected(track)) continue;\r
191         \r
192         //cut on q vectors track-by-track\r
193         if ((qPos>fTrackCuts->GetQvecMin() && qPos<fTrackCuts->GetQvecMax() && track->Eta()<0) || (qNeg>fTrackCuts->GetQvecMin() && qNeg<fTrackCuts->GetQvecMax() && track->Eta()>=0)){\r
194           \r
195           //calculate DCA for AOD track\r
196           Double_t d[2], covd[3];\r
197           AliAODTrack* track_clone=(AliAODTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters\r
198           Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);\r
199           delete track_clone;\r
200           if(!isDCA)d[0]=-999;\r
201 \r
202           fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),d[0]);  // PT histo\r
203           //Response\r
204           fHistMan->GetPIDHistogram(kHistPIDTPC)->Fill(track->GetTPCmomentum(), track->GetTPCsignal()*track->Charge()); // PID histo\r
205           fHistMan->GetPIDHistogram(kHistPIDTOF)->Fill(track->P(),track->GetTOFsignal()); // PID histo\r
206           \r
207           AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(track);\r
208           Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
209           Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)); \r
210           Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)); \r
211           Double_t nsigmaTOFkProton=0,nsigmaTOFkKaon=0,nsigmaTOFkPion=0;\r
212           if(track->Pt()>fTrackCuts->GetPtTOFMatching()){\r
213             nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
214             nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)); \r
215             nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)); \r
216             \r
217             //TOF\r
218             fHistMan->GetPtHistogram(kHistNSigProtonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
219             fHistMan->GetPtHistogram(kHistNSigKaonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));\r
220             fHistMan->GetPtHistogram(kHistNSigPionTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));\r
221             fHistMan->GetPtHistogram(kHistNSigProtonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
222             fHistMan->GetPtHistogram(kHistNSigKaonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));\r
223             fHistMan->GetPtHistogram(kHistNSigPionPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));\r
224           }\r
225           \r
226           Double_t nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);\r
227           Double_t nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);\r
228           Double_t nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);\r
229           \r
230           //TPC\r
231           fHistMan->GetPtHistogram(kHistNSigProtonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
232           fHistMan->GetPtHistogram(kHistNSigKaonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));\r
233           fHistMan->GetPtHistogram(kHistNSigPionTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));\r
234           fHistMan->GetPtHistogram(kHistNSigProtonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
235           fHistMan->GetPtHistogram(kHistNSigKaonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));\r
236           fHistMan->GetPtHistogram(kHistNSigPionPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));\r
237           //TPCTOF\r
238           fHistMan->GetPtHistogram(kHistNSigProtonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkProton);\r
239           fHistMan->GetPtHistogram(kHistNSigKaonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkKaon);\r
240           fHistMan->GetPtHistogram(kHistNSigPionTPCTOF)->Fill(track->P(),nsigmaTPCTOFkPion);\r
241          fHistMan->GetPtHistogram(kHistNSigProtonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkProton);\r
242           fHistMan->GetPtHistogram(kHistNSigKaonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkKaon);\r
243           fHistMan->GetPtHistogram(kHistNSigPionPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkPion);\r
244         \r
245           \r
246           if( ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkProton )) { \r
247             if ((nsigmaTPCTOFkKaon > fNSigmaPID) || (!CheckYCut(kKaon, track) ) ) continue;\r
248             if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaKaonPlus)->Fill(track->Pt(),d[0]); } \r
249             else { fHistMan->GetPtHistogram(kHistPtRecSigmaKaonMinus)->Fill(track->Pt(),d[0]); } \r
250           }\r
251           if( ( nsigmaTPCTOFkProton < nsigmaTPCTOFkKaon ) && ( nsigmaTPCTOFkProton < nsigmaTPCTOFkPion ) ) {\r
252             if ( nsigmaTPCTOFkProton > fNSigmaPID || (!CheckYCut(kProton, track) ) )  continue;\r
253             if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaProtonPlus)->Fill(track->Pt(),d[0]); }\r
254             else { fHistMan->GetPtHistogram(kHistPtRecSigmaProtonMinus)->Fill(track->Pt(),d[0]); }\r
255           }\r
256           if( (nsigmaTPCTOFkPion < nsigmaTPCTOFkProton ) && ( nsigmaTPCTOFkPion < nsigmaTPCTOFkKaon ) ) {\r
257             if (nsigmaTPCTOFkPion > fNSigmaPID || (!CheckYCut(kPion, track) ) ) continue;\r
258             if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaPionPlus)->Fill(track->Pt(),d[0]); }\r
259             else  { fHistMan->GetPtHistogram(kHistPtRecSigmaPionMinus)->Fill(track->Pt(),d[0]); }\r
260           }\r
261           /* MC Part */\r
262           // MF 22/02/2012\r
263           // fill mc DCA vs pt\r
264           if (arrayMC) {\r
265             AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
266             if (!partMC) { \r
267               AliError("Cannot get MC particle");\r
268               continue; }\r
269             if (CheckYCut(partMC)) {\r
270               // primaries, true pid\r
271               if ( partMC->PdgCode() == 2212) { fHistMan->GetPtHistogram(kHistPtRecTrueProtonPlus)->Fill(track->Pt(),d[0]); \r
272                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryProtonPlus)->Fill(track->Pt(),d[0]); }}\r
273               if ( partMC->PdgCode() == -2212) { fHistMan->GetPtHistogram(kHistPtRecTrueProtonMinus)->Fill(track->Pt(),d[0]); \r
274                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryProtonMinus)->Fill(track->Pt(),d[0]); }}\r
275               if ( partMC->PdgCode() == 321) { fHistMan->GetPtHistogram(kHistPtRecTrueKaonPlus)->Fill(track->Pt(),d[0]); \r
276                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryKaonPlus)->Fill(track->Pt(),d[0]); }}\r
277               if ( partMC->PdgCode() == -321) { fHistMan->GetPtHistogram(kHistPtRecTrueKaonMinus)->Fill(track->Pt(),d[0]); \r
278                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryKaonMinus)->Fill(track->Pt(),d[0]); }}\r
279               if ( partMC->PdgCode() == 211) { fHistMan->GetPtHistogram(kHistPtRecTruePionPlus)->Fill(track->Pt(),d[0]); \r
280                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionPlus)->Fill(track->Pt(),d[0]); }}\r
281               if ( partMC->PdgCode() == -211) { fHistMan->GetPtHistogram(kHistPtRecTruePionMinus)->Fill(track->Pt(),d[0]);  \r
282                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionMinus)->Fill(track->Pt(),d[0]); }}\r
283               \r
284               // primaries, sigma pid \r
285               if (partMC->IsPhysicalPrimary()) { \r
286                 if( ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkProton ) ) {\r
287                   if ( (nsigmaTPCTOFkKaon > fNSigmaPID ) || (!CheckYCut(kKaon, track) ) ) continue; \r
288                   if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryKaonPlus)->Fill(track->Pt(),d[0]); } \r
289                   else { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryKaonMinus)->Fill(track->Pt(),d[0]); } \r
290                 }\r
291                 if( ( nsigmaTPCTOFkProton < nsigmaTPCTOFkKaon ) && ( nsigmaTPCTOFkProton < nsigmaTPCTOFkPion ) ) {\r
292                   if ( (nsigmaTPCTOFkProton > fNSigmaPID ) || (!CheckYCut(kProton, track) ) ) continue;\r
293                   if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryProtonPlus)->Fill(track->Pt(),d[0]); }\r
294                   else { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryProtonMinus)->Fill(track->Pt(),d[0]); }\r
295                 }\r
296                 if( (nsigmaTPCTOFkPion < nsigmaTPCTOFkProton ) && ( nsigmaTPCTOFkPion < nsigmaTPCTOFkKaon ) ) {\r
297                   if ( ( nsigmaTPCTOFkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
298                   if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryPionPlus)->Fill(track->Pt(),d[0]); }\r
299                   else  { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryPionMinus)->Fill(track->Pt(),d[0]); }\r
300                 }\r
301               }//end if primaries\r
302               // MF 22/02/2012\r
303               // Distinguish weak decay and material\r
304               //done quickly, code can be improved\r
305               else {//to be added in a separate class\r
306                 Int_t mfl=-999,uniqueID=-999;\r
307                 Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()\r
308                 if(indexMoth>=0){//is not fake\r
309                   AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);\r
310                   Float_t codemoth = TMath::Abs(moth->GetPdgCode());\r
311                   mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));\r
312                 }\r
313                 uniqueID = partMC->GetUniqueID();\r
314                 //if(mfl==3 && uniqueID == kPDecay){//strangeness KPDecay not working on AOD, to be understood\r
315                 if(mfl==3){//strangeness\r
316                   if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) { \r
317                     if ( (nsigmaTPCkKaon > fNSigmaPID )  || (!CheckYCut(kKaon, track) ) ) continue;\r
318                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayKaonPlus)->Fill(track->Pt(),d[0]); } \r
319                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayKaonMinus)->Fill(track->Pt(),d[0]); } \r
320                   }\r
321                   if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
322                     if ( (nsigmaTPCkProton > fNSigmaPID )  || (!CheckYCut(kProton, track) ) ) continue;\r
323                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayProtonPlus)->Fill(track->Pt(),d[0]); }\r
324                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayProtonMinus)->Fill(track->Pt(),d[0]); }\r
325                   }\r
326                   if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
327                     if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
328                     if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayPionPlus)->Fill(track->Pt(),d[0]); }\r
329                     else  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayPionMinus)->Fill(track->Pt(),d[0]); }\r
330                   }\r
331                 }//end if strangeness\r
332                 else{//material\r
333                   if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) { \r
334                     if ( (nsigmaTPCkKaon > fNSigmaPID )  || (!CheckYCut(kKaon, track) ) ) continue;\r
335                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialKaonPlus)->Fill(track->Pt(),d[0]); } \r
336                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialKaonMinus)->Fill(track->Pt(),d[0]); } \r
337                   }\r
338                   if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
339                     if ( (nsigmaTPCkProton > fNSigmaPID )  || (!CheckYCut(kProton, track) ) ) continue;\r
340                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialProtonPlus)->Fill(track->Pt(),d[0]); }\r
341                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialProtonMinus)->Fill(track->Pt(),d[0]); }\r
342                   }\r
343                   if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
344                     if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
345                     if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialPionPlus)->Fill(track->Pt(),d[0]); }\r
346                     else  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialPionMinus)->Fill(track->Pt(),d[0]); }\r
347                   }\r
348                 }//end if material\r
349               }//end else (IsPhysPrim)\r
350             }//end if(CheckY)\r
351           }//end if(arrayMC)\r
352         }//end if on q vector (track)\r
353       } // end loop on tracks\r
354   }//end if on q vector (event)\r
355   PostData(1, fHistMan);\r
356   PostData(2, fEventCuts);\r
357   PostData(3, fTrackCuts);\r
358 }\r
359 \r
360 //_________________________________________________________________\r
361 void   AliAnalysisTaskSpectraAOD::Terminate(Option_t *)\r
362 {\r
363   // Terminate\r
364 }\r