cff85b53c0016f28f4e62fe06b0bbdcf6a563875
[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 == kProton) { y = track->Y(9.38271999999999995e-01); }\r
74   if ( species == kKaon ) { y = track->Y(4.93676999999999977e-01); }\r
75   if ( species == kPion)  { 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(TMath::Abs(partMC->Eta()) > fTrackCuts->GetEta()) continue;\r
169             \r
170             fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
171             // check for true PID + and fill P_t histos \r
172             //if (partMC->IsPhysicalPrimary() && CheckYCut(partMC) ) {// only primary vertices and y cut satisfied\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     \r
185     \r
186     //main loop on tracks\r
187     for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)\r
188       {\r
189         \r
190         AliAODTrack* track = fAOD->GetTrack(iTracks);\r
191         if (!fTrackCuts->IsSelected(track)) continue;\r
192         \r
193         //cut on q vectors track-by-track\r
194         if ((qPos>fTrackCuts->GetQvecMin() && qPos<fTrackCuts->GetQvecMax() && track->Eta()<0) || (qNeg>fTrackCuts->GetQvecMin() && qNeg<fTrackCuts->GetQvecMax() && track->Eta()>=0)){\r
195           \r
196           //calculate DCA for AOD track\r
197           Double_t d[2], covd[3];\r
198           AliAODTrack* track_clone=(AliAODTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters\r
199           Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);\r
200           delete track_clone;\r
201           if(!isDCA)d[0]=-999;\r
202           \r
203           fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),d[0]);  // PT histo\r
204           //Response\r
205           fHistMan->GetPIDHistogram(kHistPIDTPC)->Fill(track->GetTPCmomentum(), track->GetTPCsignal()*track->Charge()); // 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             fHistMan->GetPIDHistogram(kHistPIDTOF)->Fill(track->P(),(track->GetTOFsignal()/100)*track->Charge()); // PID histo\r
214             nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
215             nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)); \r
216             nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)); \r
217             \r
218             //TOF\r
219             fHistMan->GetPtHistogram(kHistNSigProtonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
220             fHistMan->GetPtHistogram(kHistNSigKaonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));\r
221             fHistMan->GetPtHistogram(kHistNSigPionTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));\r
222             fHistMan->GetPtHistogram(kHistNSigProtonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));\r
223             fHistMan->GetPtHistogram(kHistNSigKaonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));\r
224             fHistMan->GetPtHistogram(kHistNSigPionPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));\r
225           }\r
226           \r
227           Double_t nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);\r
228           Double_t nsigmaTPCTOFkKaon = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);\r
229           Double_t nsigmaTPCTOFkPion = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);\r
230           \r
231           //TPC\r
232           fHistMan->GetPtHistogram(kHistNSigProtonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
233           fHistMan->GetPtHistogram(kHistNSigKaonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));\r
234           fHistMan->GetPtHistogram(kHistNSigPionTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));\r
235           fHistMan->GetPtHistogram(kHistNSigProtonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
236           fHistMan->GetPtHistogram(kHistNSigKaonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon));\r
237           fHistMan->GetPtHistogram(kHistNSigPionPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion));\r
238           //TPCTOF\r
239           fHistMan->GetPtHistogram(kHistNSigProtonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkProton);\r
240           fHistMan->GetPtHistogram(kHistNSigKaonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkKaon);\r
241           fHistMan->GetPtHistogram(kHistNSigPionTPCTOF)->Fill(track->P(),nsigmaTPCTOFkPion);\r
242           fHistMan->GetPtHistogram(kHistNSigProtonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkProton);\r
243           fHistMan->GetPtHistogram(kHistNSigKaonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkKaon);\r
244           fHistMan->GetPtHistogram(kHistNSigPionPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkPion);\r
245         \r
246           \r
247           if( ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkProton )) { \r
248             if ((nsigmaTPCTOFkKaon > fNSigmaPID) || (!CheckYCut(kKaon, track) ) ) continue;\r
249             if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaKaonPlus)->Fill(track->Pt(),d[0]); } \r
250             else { fHistMan->GetPtHistogram(kHistPtRecSigmaKaonMinus)->Fill(track->Pt(),d[0]); } \r
251           }\r
252           if( ( nsigmaTPCTOFkProton < nsigmaTPCTOFkKaon ) && ( nsigmaTPCTOFkProton < nsigmaTPCTOFkPion ) ) {\r
253             if ( nsigmaTPCTOFkProton > fNSigmaPID || (!CheckYCut(kProton, track) ) )  continue;\r
254             if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaProtonPlus)->Fill(track->Pt(),d[0]); }\r
255             else { fHistMan->GetPtHistogram(kHistPtRecSigmaProtonMinus)->Fill(track->Pt(),d[0]); }\r
256           }\r
257           if( (nsigmaTPCTOFkPion < nsigmaTPCTOFkProton ) && ( nsigmaTPCTOFkPion < nsigmaTPCTOFkKaon ) ) {\r
258             if (nsigmaTPCTOFkPion > fNSigmaPID || (!CheckYCut(kPion, track) ) ) continue;\r
259             if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaPionPlus)->Fill(track->Pt(),d[0]); }\r
260             else  { fHistMan->GetPtHistogram(kHistPtRecSigmaPionMinus)->Fill(track->Pt(),d[0]); }\r
261           }\r
262           /* MC Part */\r
263           // MF 22/02/2012\r
264           // fill mc DCA vs pt\r
265           if (arrayMC) {\r
266             AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
267             if (!partMC) { \r
268               AliError("Cannot get MC particle");\r
269               continue; }\r
270             if (partMC->IsPhysicalPrimary())fHistMan->GetPtHistogram(kHistPtRecPrimary)->Fill(track->Pt(),d[0]);  // PT histo\r
271             if (CheckYCut(partMC)){\r
272               // primaries, true pid\r
273               //25th Apr - nsigma cut added in addition to the PDG code\r
274               if ( partMC->PdgCode() == 2212 && nsigmaTPCTOFkProton < fNSigmaPID) { fHistMan->GetPtHistogram(kHistPtRecTrueProtonPlus)->Fill(track->Pt(),d[0]); \r
275                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryProtonPlus)->Fill(track->Pt(),d[0]); }}\r
276               if ( partMC->PdgCode() == -2212 && nsigmaTPCTOFkProton < fNSigmaPID) { fHistMan->GetPtHistogram(kHistPtRecTrueProtonMinus)->Fill(track->Pt(),d[0]); \r
277                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryProtonMinus)->Fill(track->Pt(),d[0]); }}\r
278               if ( partMC->PdgCode() == 321 && nsigmaTPCTOFkKaon < fNSigmaPID) { fHistMan->GetPtHistogram(kHistPtRecTrueKaonPlus)->Fill(track->Pt(),d[0]); \r
279                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryKaonPlus)->Fill(track->Pt(),d[0]); }}\r
280               if ( partMC->PdgCode() == -321 && nsigmaTPCTOFkKaon < fNSigmaPID) { fHistMan->GetPtHistogram(kHistPtRecTrueKaonMinus)->Fill(track->Pt(),d[0]); \r
281                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryKaonMinus)->Fill(track->Pt(),d[0]); }}\r
282               if ( partMC->PdgCode() == 211 && nsigmaTPCTOFkPion < fNSigmaPID) { fHistMan->GetPtHistogram(kHistPtRecTruePionPlus)->Fill(track->Pt(),d[0]); \r
283                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionPlus)->Fill(track->Pt(),d[0]); }}\r
284               if ( partMC->PdgCode() == -211 && nsigmaTPCTOFkPion < fNSigmaPID) { fHistMan->GetPtHistogram(kHistPtRecTruePionMinus)->Fill(track->Pt(),d[0]);  \r
285                 if (partMC->IsPhysicalPrimary()) {fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionMinus)->Fill(track->Pt(),d[0]); }}\r
286               //25th Apr - Muons are added to Pions\r
287               if ( partMC->PdgCode() == 13 && nsigmaTPCTOFkPion < fNSigmaPID) { \r
288                 fHistMan->GetPtHistogram(kHistPtRecTruePionPlus)->Fill(track->Pt(),d[0]); \r
289                 fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionPlus)->Fill(track->Pt(),d[0]);///////////////////FIXME \r
290                 fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),d[0]); \r
291                 if (partMC->IsPhysicalPrimary()) {\r
292                   fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),d[0]); \r
293                 }}\r
294               if ( partMC->PdgCode() == -13 && nsigmaTPCTOFkPion < fNSigmaPID) { \r
295                 fHistMan->GetPtHistogram(kHistPtRecTruePionMinus)->Fill(track->Pt(),d[0]); \r
296                 fHistMan->GetPtHistogram(kHistPtRecTruePrimaryPionMinus)->Fill(track->Pt(),d[0]);//////////////FIXME \r
297                 fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),d[0]); \r
298                 if (partMC->IsPhysicalPrimary()) {\r
299                   fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),d[0]); \r
300                 }}\r
301               \r
302               // primaries, sigma pid \r
303               if (partMC->IsPhysicalPrimary()) { \r
304                 if( ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkPion ) && ( nsigmaTPCTOFkKaon < nsigmaTPCTOFkProton ) ) {\r
305                   if ( (nsigmaTPCTOFkKaon > fNSigmaPID ) || (!CheckYCut(kKaon, track) ) ) continue; \r
306                   if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryKaonPlus)->Fill(track->Pt(),d[0]); } \r
307                   else { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryKaonMinus)->Fill(track->Pt(),d[0]); } \r
308                 }\r
309                 if( ( nsigmaTPCTOFkProton < nsigmaTPCTOFkKaon ) && ( nsigmaTPCTOFkProton < nsigmaTPCTOFkPion ) ) {\r
310                   if ( (nsigmaTPCTOFkProton > fNSigmaPID ) || (!CheckYCut(kProton, track) ) ) continue;\r
311                   if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryProtonPlus)->Fill(track->Pt(),d[0]); }\r
312                   else { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryProtonMinus)->Fill(track->Pt(),d[0]); }\r
313                 }\r
314                 if( (nsigmaTPCTOFkPion < nsigmaTPCTOFkProton ) && ( nsigmaTPCTOFkPion < nsigmaTPCTOFkKaon ) ) {\r
315                   if ( ( nsigmaTPCTOFkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
316                   if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryPionPlus)->Fill(track->Pt(),d[0]); }\r
317                   else  { fHistMan->GetPtHistogram(kHistPtRecSigmaPrimaryPionMinus)->Fill(track->Pt(),d[0]); }\r
318                 }\r
319               }//end if primaries\r
320               // MF 22/02/2012\r
321               // Distinguish weak decay and material\r
322               //done quickly, code can be improved\r
323               else {//to be added in a separate class\r
324                 Int_t mfl=-999,codemoth=-999;\r
325                 Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()\r
326                 if(indexMoth>=0){//is not fake\r
327                   AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);\r
328                   codemoth = TMath::Abs(moth->GetPdgCode());\r
329                   mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));\r
330                 }\r
331                 //UInt_t flag; \r
332                 //flag=partMC->GetStatus();\r
333                 //Printf("FLAG: %d",flag);\r
334                 //if(mfl==3 && (flag&AliAODMCParticle::kPDecay)!=0){//strangeness KPDecay not working on AOD, to be understood\r
335                 //Printf("\n\n\n STRANGENESS!!!!!");\r
336                 //if(codemoth!=-999)Printf("mfl:%d    codemoth%d",mfl,codemoth);\r
337                 //if(codemoth==310 || codemoth==130 || codemoth==321 || codemoth==3122 || codemoth==3112 ||\r
338                 //codemoth==3222 || codemoth==3312 || codemoth==3322 || codemoth==3334){//K0_S, K0_L, K^+-,lambda, sigma0,sigma+,xi-,xi0, omega\r
339                 if(mfl==3){//strangeness\r
340                   if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) { \r
341                     if ( (nsigmaTPCkKaon > fNSigmaPID )  || (!CheckYCut(kKaon, track) ) ) continue;\r
342                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayKaonPlus)->Fill(track->Pt(),d[0]); } \r
343                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayKaonMinus)->Fill(track->Pt(),d[0]); } \r
344                   }\r
345                   if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
346                     if ( (nsigmaTPCkProton > fNSigmaPID )  || (!CheckYCut(kProton, track) ) ) continue;\r
347                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayProtonPlus)->Fill(track->Pt(),d[0]); }\r
348                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayProtonMinus)->Fill(track->Pt(),d[0]); }\r
349                   }\r
350                   if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
351                     if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
352                     if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayPionPlus)->Fill(track->Pt(),d[0]); }\r
353                     else  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryWeakDecayPionMinus)->Fill(track->Pt(),d[0]); }\r
354                   }\r
355                 }//end if strangeness\r
356                 else{//material\r
357                   if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) { \r
358                     if ( (nsigmaTPCkKaon > fNSigmaPID )  || (!CheckYCut(kKaon, track) ) ) continue;\r
359                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialKaonPlus)->Fill(track->Pt(),d[0]); } \r
360                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialKaonMinus)->Fill(track->Pt(),d[0]); } \r
361                   }\r
362                   if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
363                     if ( (nsigmaTPCkProton > fNSigmaPID )  || (!CheckYCut(kProton, track) ) ) continue;\r
364                     if ( track->Charge() > 0 ) { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialProtonPlus)->Fill(track->Pt(),d[0]); }\r
365                     else { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialProtonMinus)->Fill(track->Pt(),d[0]); }\r
366                   }\r
367                   if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
368                     if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
369                     if ( track->Charge() > 0 )  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialPionPlus)->Fill(track->Pt(),d[0]); }\r
370                     else  { fHistMan->GetPtHistogram(kHistPtRecSigmaSecondaryMaterialPionMinus)->Fill(track->Pt(),d[0]); }\r
371                   }\r
372                 }//end if material\r
373               }//end else (IsPhysPrim)\r
374             }//end if(CheckY)\r
375           }//end if(arrayMC)\r
376         }//end if on q vector (track)\r
377       } // end loop on tracks\r
378   }//end if on q vector (event)\r
379   PostData(1, fHistMan);\r
380   PostData(2, fEventCuts);\r
381   PostData(3, fTrackCuts);\r
382 }\r
383 \r
384 //_________________________________________________________________\r
385 void   AliAnalysisTaskSpectraAOD::Terminate(Option_t *)\r
386 {\r
387   // Terminate\r
388 }\r