Implemented handy getters + renamed some of the enums to avoid conflicts with root...
[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 "AliVParticle.h"\r
31 #include "AliAODEvent.h"\r
32 #include "AliAODInputHandler.h"\r
33 #include "AliAnalysisTaskSpectraAOD.h"\r
34 #include "AliAnalysisTaskESDfilter.h"\r
35 #include "AliAnalysisDataContainer.h"\r
36 #include "AliSpectraAODHistoManager.h"\r
37 #include "AliSpectraAODTrackCuts.h"\r
38 #include "AliSpectraAODEventCuts.h"\r
39 #include "AliCentrality.h"\r
40 #include "TProof.h"\r
41 #include "AliPID.h"\r
42 #include "AliVEvent.h"\r
43 #include "AliPIDResponse.h"\r
44 #include "AliStack.h"\r
45 #include <iostream>\r
46 \r
47 \r
48 \r
49 \r
50 using namespace AliSpectraNameSpace;\r
51 using namespace std;\r
52 \r
53 ClassImp(AliAnalysisTaskSpectraAOD) // EX1 This stuff tells root to implement the streamer, inspector methods etc (we discussed about it today)\r
54 \r
55 //________________________________________________________________________\r
56 AliAnalysisTaskSpectraAOD::AliAnalysisTaskSpectraAOD(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0), fIsMC(0), fPIDResponse(0), fNSigmaPID(0), fYCut(0)\r
57 {\r
58   // Default constructor\r
59   fNSigmaPID = 3.;\r
60   fYCut = .5;\r
61   DefineInput(0, TChain::Class());\r
62   DefineOutput(1, AliSpectraAODHistoManager::Class());\r
63   DefineOutput(2, AliSpectraAODEventCuts::Class());\r
64   DefineOutput(3, AliSpectraAODTrackCuts::Class());\r
65 \r
66 }\r
67 //________________________________________________________________________\r
68 Bool_t AliAnalysisTaskSpectraAOD::CheckYCut(AODParticleSpecies_t species, AliAODTrack* track) const\r
69 {\r
70   // check if the rapidity is within the set range\r
71   // 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
72   Double_t y;\r
73   if (species == kSpProton) { y = track->Y(9.38271999999999995e-01); }\r
74   if ( species == kSpKaon ) { y = track->Y(4.93676999999999977e-01); }\r
75   if ( species == kSpPion)  { y = track->Y(1.39570000000000000e-01); }\r
76   if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;\r
77   return kTRUE;\r
78 }\r
79 //____________________________________________________________________________\r
80 Bool_t AliAnalysisTaskSpectraAOD::CheckYCut(AliAODMCParticle* particle) const\r
81 {\r
82   // check if the rapidity is within the set range\r
83   Double_t y = particle->Y();\r
84   if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;\r
85   return kTRUE;\r
86 }\r
87 //________________________________________________________________________\r
88 void AliAnalysisTaskSpectraAOD::UserCreateOutputObjects()\r
89 {\r
90   // create output objects\r
91   fHistMan = new AliSpectraAODHistoManager("SpectraHistos");\r
92 \r
93   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
94   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
95 \r
96   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();\r
97   AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());\r
98   fPIDResponse = inputHandler->GetPIDResponse();\r
99 \r
100   PostData(1, fHistMan);\r
101   PostData(2, fEventCuts);\r
102   PostData(3, fTrackCuts);\r
103 \r
104 }\r
105 //________________________________________________________________________\r
106 void AliAnalysisTaskSpectraAOD::UserExec(Option_t *)\r
107 {\r
108   // main event loop\r
109   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);\r
110   if (strcmp(fAOD->ClassName(), "AliAODEvent"))\r
111     {\r
112       AliFatal("Not processing AODs");\r
113     }\r
114   \r
115   if(!fPIDResponse) {\r
116     AliError("Cannot get pid response");\r
117     return;\r
118   }\r
119   \r
120   //Selection on QVector, before ANY other selection on the event\r
121   //Spectra MUST be normalized wrt events AFTER the selection on Qvector\r
122   // Can we include this in fEventCuts\r
123   Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;\r
124   Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;\r
125   Int_t multPos = 0;\r
126   Int_t multNeg = 0;\r
127   for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {\r
128     AliAODTrack* aodTrack = fAOD->GetTrack(iT);\r
129     if (!fTrackCuts->IsSelected(aodTrack)) continue;\r
130     if (aodTrack->Eta() >= 0){\r
131       multPos++;\r
132       Qx2EtaPos += TMath::Cos(2*aodTrack->Phi()); \r
133       Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());\r
134     } else {\r
135       multNeg++;\r
136       Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi()); \r
137       Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());\r
138     }\r
139   } \r
140   Double_t qPos = TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);\r
141   Double_t qNeg = TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);\r
142   if((qPos>fTrackCuts->GetQvecMin() && qPos<fTrackCuts->GetQvecMax()) || (qNeg>fTrackCuts->GetQvecMin() && qNeg<fTrackCuts->GetQvecMax())){\r
143     \r
144     //check on centrality distribution\r
145     fHistMan->GetPtHistogram("CentCheck")->Fill(fAOD->GetCentrality()->GetCentralityPercentile("V0M"),fAOD->GetHeader()->GetCentralityP()->GetCentralityPercentileUnchecked("V0M"));\r
146     \r
147     if (!fEventCuts->IsSelected(fAOD)) return;\r
148     \r
149     //fill q distributions vs centrality, after all event selection\r
150     fHistMan->GetqVecHistogram(kHistqVecPos)->Fill(qPos,fAOD->GetCentrality()->GetCentralityPercentile("V0M"));  // qVector distribution\r
151     fHistMan->GetqVecHistogram(kHistqVecNeg)->Fill(qNeg,fAOD->GetCentrality()->GetCentralityPercentile("V0M"));  // qVector distribution\r
152     \r
153     //AliCentrality fAliCentral*;\r
154     //   if ((fAOD->GetCentrality())->GetCentralityPercentile("V0M") > 5.) return;\r
155     \r
156     // First do MC to fill up the MC particle array, such that we can use it later\r
157     TClonesArray *arrayMC = 0;\r
158     if (fIsMC)\r
159       {\r
160         arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
161         if (!arrayMC) {\r
162           AliFatal("Error: MC particles branch not found!\n");\r
163         }\r
164         Int_t nMC = arrayMC->GetEntries();\r
165         for (Int_t iMC = 0; iMC < nMC; iMC++)\r
166           {\r
167             AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
168             if(!partMC->Charge()) continue;//Skip neutrals\r
169             if(TMath::Abs(partMC->Eta()) > fTrackCuts->GetEta()) continue;\r
170             \r
171             fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
172             // check for true PID + and fill P_t histos \r
173             if (CheckYCut(partMC) ){// only primary vertices and y cut satisfied\r
174               if ( partMC->PdgCode() == 2212) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryProtonPlus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
175               if ( partMC->PdgCode() == -2212) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryProtonMinus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
176               if ( partMC->PdgCode() == 321)  { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryKaonPlus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
177               if ( partMC->PdgCode() == -321) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryKaonMinus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
178               if ( partMC->PdgCode() == 211)  { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryPionPlus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
179               if ( partMC->PdgCode() == -211) { fHistMan->GetPtHistogram(kHistPtGenTruePrimaryPionMinus)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } \r
180             }\r
181           }\r
182       }\r
183     \r
184     //main loop on tracks\r
185     for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)\r
186       {\r
187         AliAODTrack* track = fAOD->GetTrack(iTracks);\r
188         if (!fTrackCuts->IsSelected(track)) continue;\r
189         \r
190         //cut on q vectors track-by-track\r
191         if ((qPos>fTrackCuts->GetQvecMin() && qPos<fTrackCuts->GetQvecMax() && track->Eta()<0) || (qNeg>fTrackCuts->GetQvecMin() && qNeg<fTrackCuts->GetQvecMax() && track->Eta()>=0)){\r
192           \r
193           //calculate DCA for AOD track\r
194           Double_t d[2], covd[3];\r
195           AliAODTrack* track_clone=(AliAODTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters\r
196           Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);\r
197           delete track_clone;\r
198           if(!isDCA)d[0]=-999;\r
199           \r
200           fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),d[0]);  // PT histo\r
201           //Response\r
202           fHistMan->GetPIDHistogram(kHistPIDTPC)->Fill(track->GetTPCmomentum(), track->GetTPCsignal()*track->Charge()); // PID histo\r
203           \r
204           AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(track);\r
205           Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
206           Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)); \r
207           Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)); \r
208           Double_t nsigmaTOFkProton=0,nsigmaTOFkKaon=0,nsigmaTOFkPion=0;\r
209           if(track->Pt()>fTrackCuts->GetPtTOFMatching()){\r
210             fHistMan->GetPIDHistogram(kHistPIDTOF)->Fill(track->P(),(track->GetTOFsignal()/100)*track->Charge()); // PID histo\r
211             nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
212             nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)); \r
213             nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)); \r
214             \r
215             //TOF\r
216             fHistMan->GetPtHistogram(kHistNSigProtonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
217             fHistMan->GetPtHistogram(kHistNSigKaonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));\r
218             fHistMan->GetPtHistogram(kHistNSigPionTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));\r
219             fHistMan->GetPtHistogram(kHistNSigProtonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
220             fHistMan->GetPtHistogram(kHistNSigKaonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));\r
221             fHistMan->GetPtHistogram(kHistNSigPionPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));\r
222           }\r
223           \r
224           Double_t nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);\r
225           Double_t nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);\r
226           Double_t nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);\r
227           \r
228           //TPC\r
229           fHistMan->GetPtHistogram(kHistNSigProtonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
230           fHistMan->GetPtHistogram(kHistNSigKaonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));\r
231           fHistMan->GetPtHistogram(kHistNSigPionTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));\r
232           fHistMan->GetPtHistogram(kHistNSigProtonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
233           fHistMan->GetPtHistogram(kHistNSigKaonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));\r
234           fHistMan->GetPtHistogram(kHistNSigPionPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));\r
235           //TPCTOF\r
236           fHistMan->GetPtHistogram(kHistNSigProtonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkProton);\r
237           fHistMan->GetPtHistogram(kHistNSigKaonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkKaon);\r
238           fHistMan->GetPtHistogram(kHistNSigPionTPCTOF)->Fill(track->P(),nsigmaTPCTOFkPion);\r
239           fHistMan->GetPtHistogram(kHistNSigProtonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkProton);\r
240           fHistMan->GetPtHistogram(kHistNSigKaonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkKaon);\r
241           fHistMan->GetPtHistogram(kHistNSigPionPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkPion);\r
242         \r
243           \r
244           if( ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkProton )) { \r
245             if ((nsigmaTPCTOFkKaon > fNSigmaPID) || (!CheckYCut(kSpKaon, track) ) ) continue;\r
246             if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaKaonPlus)->Fill(track->Pt(),d[0]); } \r
247             else { fHistMan->GetPtHistogram(kHistPtRecSigmaKaonMinus)->Fill(track->Pt(),d[0]); } \r
248           }\r
249           if( ( nsigmaTPCTOFkProton < nsigmaTPCTOFkKaon ) && ( nsigmaTPCTOFkProton < nsigmaTPCTOFkPion ) ) {\r
250             if ( nsigmaTPCTOFkProton > fNSigmaPID || (!CheckYCut(kSpProton, track) ) )  continue;\r
251             if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaProtonPlus)->Fill(track->Pt(),d[0]); }\r
252             else { fHistMan->GetPtHistogram(kHistPtRecSigmaProtonMinus)->Fill(track->Pt(),d[0]); }\r
253           }\r
254           if( (nsigmaTPCTOFkPion < nsigmaTPCTOFkProton ) && ( nsigmaTPCTOFkPion < nsigmaTPCTOFkKaon ) ) {\r
255             if (nsigmaTPCTOFkPion > fNSigmaPID || (!CheckYCut(kSpPion, track) ) ) continue;\r
256             if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaPionPlus)->Fill(track->Pt(),d[0]); }\r
257             else  { fHistMan->GetPtHistogram(kHistPtRecSigmaPionMinus)->Fill(track->Pt(),d[0]); }\r
258           }\r
259           /* MC Part */\r
260           // MF 22/02/2012\r
261           // fill mc DCA vs pt\r
262           if (arrayMC) {\r
263             AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
264             if (!partMC) { \r
265               AliError("Cannot get MC particle");\r
266               continue; }\r
267             if (partMC->IsPhysicalPrimary())fHistMan->GetPtHistogram(kHistPtRecPrimary)->Fill(track->Pt(),d[0]);  // PT histo\r
268             if (CheckYCut(partMC)){\r
269               // primaries, true pid\r
270               //25th Apr - nsigma cut added in addition to the PDG code\r
271               if ( partMC->PdgCode() == 2212 && nsigmaTPCTOFkProton < fNSigmaPID) { 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 && nsigmaTPCTOFkProton < fNSigmaPID) { 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 && nsigmaTPCTOFkKaon < fNSigmaPID) { 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 && nsigmaTPCTOFkKaon < fNSigmaPID) { 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 && nsigmaTPCTOFkPion < fNSigmaPID) { 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 && nsigmaTPCTOFkPion < fNSigmaPID) { fHistMan->GetPtHistogram(kHistPtRecTruePionMinus)->Fill(track->Pt(),d[0]);  \r
282                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionMinus)->Fill(track->Pt(),d[0]); }}\r
283               //25th Apr - Muons are added to Pions\r
284               if ( partMC->PdgCode() == 13 && nsigmaTPCTOFkPion < fNSigmaPID) { \r
285                 fHistMan->GetPtHistogram(kHistPtRecTruePionPlus)->Fill(track->Pt(),d[0]); \r
286                 fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionPlus)->Fill(track->Pt(),d[0]);///////////////////FIXME \r
287                 fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),d[0]); \r
288                 if (partMC->IsPhysicalPrimary()) {\r
289                   fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),d[0]); \r
290                 }}\r
291               if ( partMC->PdgCode() == -13 && nsigmaTPCTOFkPion < fNSigmaPID) { \r
292                 fHistMan->GetPtHistogram(kHistPtRecTruePionMinus)->Fill(track->Pt(),d[0]); \r
293                 fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionMinus)->Fill(track->Pt(),d[0]);//////////////FIXME \r
294                 fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),d[0]); \r
295                 if (partMC->IsPhysicalPrimary()) {\r
296                   fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),d[0]); \r
297                 }}\r
298               \r
299               // primaries, sigma pid \r
300               if (partMC->IsPhysicalPrimary()) { \r
301                 if( ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkProton ) ) {\r
302                   if ( (nsigmaTPCTOFkKaon > fNSigmaPID ) || (!CheckYCut(kSpKaon, track) ) ) continue; \r
303                   if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryKaonPlus)->Fill(track->Pt(),d[0]); } \r
304                   else { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryKaonMinus)->Fill(track->Pt(),d[0]); } \r
305                 }\r
306                 if( ( nsigmaTPCTOFkProton < nsigmaTPCTOFkKaon ) && ( nsigmaTPCTOFkProton < nsigmaTPCTOFkPion ) ) {\r
307                   if ( (nsigmaTPCTOFkProton > fNSigmaPID ) || (!CheckYCut(kSpProton, track) ) ) continue;\r
308                   if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryProtonPlus)->Fill(track->Pt(),d[0]); }\r
309                   else { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryProtonMinus)->Fill(track->Pt(),d[0]); }\r
310                 }\r
311                 if( (nsigmaTPCTOFkPion < nsigmaTPCTOFkProton ) && ( nsigmaTPCTOFkPion < nsigmaTPCTOFkKaon ) ) {\r
312                   if ( ( nsigmaTPCTOFkPion > fNSigmaPID )  || (!CheckYCut(kSpPion, track) ) ) continue;\r
313                   if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryPionPlus)->Fill(track->Pt(),d[0]); }\r
314                   else  { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryPionMinus)->Fill(track->Pt(),d[0]); }\r
315                 }\r
316               }//end if primaries\r
317               // MF 22/02/2012\r
318               // Distinguish weak decay and material\r
319               //done quickly, code can be improved\r
320               else {//to be added in a separate class\r
321                 Int_t mfl=-999,codemoth=-999;\r
322                 Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()\r
323                 if(indexMoth>=0){//is not fake\r
324                   AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);\r
325                   codemoth = TMath::Abs(moth->GetPdgCode());\r
326                   mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));\r
327                 }\r
328                 //UInt_t flag; \r
329                 //flag=partMC->GetStatus();\r
330                 //Printf("FLAG: %d",flag);\r
331                 //if(mfl==3 && (flag&AliAODMCParticle::kPDecay)!=0){//strangeness KPDecay not working on AOD, to be understood\r
332                 //Printf("\n\n\n STRANGENESS!!!!!");\r
333                 //if(codemoth!=-999)Printf("mfl:%d    codemoth%d",mfl,codemoth);\r
334                 //if(codemoth==310 || codemoth==130 || codemoth==321 || codemoth==3122 || codemoth==3112 ||\r
335                 //codemoth==3222 || codemoth==3312 || codemoth==3322 || codemoth==3334){//K0_S, K0_L, K^+-,lambda, sigma0,sigma+,xi-,xi0, omega\r
336                 if(mfl==3){//strangeness\r
337                   if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) { \r
338                     if ( (nsigmaTPCkKaon > fNSigmaPID )  || (!CheckYCut(kSpKaon, track) ) ) continue;\r
339                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayKaonPlus)->Fill(track->Pt(),d[0]); } \r
340                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayKaonMinus)->Fill(track->Pt(),d[0]); } \r
341                   }\r
342                   if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
343                     if ( (nsigmaTPCkProton > fNSigmaPID )  || (!CheckYCut(kSpProton, track) ) ) continue;\r
344                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayProtonPlus)->Fill(track->Pt(),d[0]); }\r
345                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayProtonMinus)->Fill(track->Pt(),d[0]); }\r
346                   }\r
347                   if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
348                     if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kSpPion, track) ) ) continue;\r
349                     if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayPionPlus)->Fill(track->Pt(),d[0]); }\r
350                     else  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayPionMinus)->Fill(track->Pt(),d[0]); }\r
351                   }\r
352                 }//end if strangeness\r
353                 else{//material\r
354                   if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) { \r
355                     if ( (nsigmaTPCkKaon > fNSigmaPID )  || (!CheckYCut(kSpKaon, track) ) ) continue;\r
356                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialKaonPlus)->Fill(track->Pt(),d[0]); } \r
357                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialKaonMinus)->Fill(track->Pt(),d[0]); } \r
358                   }\r
359                   if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
360                     if ( (nsigmaTPCkProton > fNSigmaPID )  || (!CheckYCut(kSpProton, track) ) ) continue;\r
361                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialProtonPlus)->Fill(track->Pt(),d[0]); }\r
362                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialProtonMinus)->Fill(track->Pt(),d[0]); }\r
363                   }\r
364                   if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
365                     if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kSpPion, track) ) ) continue;\r
366                     if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialPionPlus)->Fill(track->Pt(),d[0]); }\r
367                     else  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialPionMinus)->Fill(track->Pt(),d[0]); }\r
368                   }\r
369                 }//end if material\r
370               }//end else (IsPhysPrim)\r
371             }//end if(CheckY)\r
372           }//end if(arrayMC)\r
373         }//end if on q vector (track)\r
374       } // end loop on tracks\r
375   }//end if on q vector (event)\r
376   PostData(1, fHistMan);\r
377   PostData(2, fEventCuts);\r
378   PostData(3, fTrackCuts);\r
379 }\r
380 \r
381 //_________________________________________________________________\r
382 void   AliAnalysisTaskSpectraAOD::Terminate(Option_t *)\r
383 {\r
384   // Terminate\r
385 }\r