- processed events + QA plots on vtx and dndeta added in Event Cuts - small bugfix...
[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 "AliSpectraAODPID.h"\r
46 #include <iostream>\r
47 \r
48 \r
49 \r
50 \r
51 using namespace AliSpectraNameSpace;\r
52 using namespace std;\r
53 \r
54 ClassImp(AliAnalysisTaskSpectraAOD) // EX1 This stuff tells root to implement the streamer, inspector methods etc (we discussed about it today)\r
55 \r
56 //________________________________________________________________________\r
57 AliAnalysisTaskSpectraAOD::AliAnalysisTaskSpectraAOD(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0),  fPID(0), fIsMC(0)\r
58 {\r
59   // Default constructor\r
60 \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   DefineOutput(4, AliSpectraAODPID::Class());\r
66 \r
67 }\r
68 //________________________________________________________________________\r
69 //________________________________________________________________________\r
70 void AliAnalysisTaskSpectraAOD::UserCreateOutputObjects()\r
71 {\r
72   // create output objects\r
73   fHistMan = new AliSpectraAODHistoManager("SpectraHistos");\r
74 \r
75   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
76   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
77   if (!fPID)       AliFatal("PID object should be set in the steering macro");\r
78 \r
79   PostData(1, fHistMan  );\r
80   PostData(2, fEventCuts);\r
81   PostData(3, fTrackCuts);\r
82   PostData(4, fPID      );\r
83 \r
84 }\r
85 //________________________________________________________________________\r
86 void AliAnalysisTaskSpectraAOD::UserExec(Option_t *)\r
87 {\r
88   // main event loop\r
89   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);\r
90   if (strcmp(fAOD->ClassName(), "AliAODEvent"))\r
91     {\r
92       AliFatal("Not processing AODs");\r
93     }\r
94     \r
95   //check on centrality distribution\r
96   fHistMan->GetPtHistogram("CentCheck")->Fill(fAOD->GetCentrality()->GetCentralityPercentile("V0M"),fAOD->GetHeader()->GetCentralityP()->GetCentralityPercentileUnchecked("V0M"));\r
97   \r
98   if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
99   \r
100   //AliCentrality fAliCentral*;\r
101   //   if ((fAOD->GetCentrality())->GetCentralityPercentile("V0M") > 5.) return;\r
102   \r
103   // First do MC to fill up the MC particle array, such that we can use it later\r
104   TClonesArray *arrayMC = 0;\r
105   if (fIsMC)\r
106     {\r
107       arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
108       if (!arrayMC) {\r
109         AliFatal("Error: MC particles branch not found!\n");\r
110       }\r
111       Int_t nMC = arrayMC->GetEntries();\r
112       for (Int_t iMC = 0; iMC < nMC; iMC++)\r
113         {\r
114           AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
115           if(!partMC->Charge()) continue;//Skip neutrals\r
116           //if(TMath::Abs(partMC->Eta()) > fTrackCuts->GetEta()) continue; FIXME eta cut should not be there!!!! we want the generated in y<0.5\r
117           fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
118           \r
119           if(TMath::Abs(partMC->Y())   > fTrackCuts->GetY()  ) continue;            \r
120           // check for true PID + and fill P_t histos \r
121           Int_t charge = partMC->Charge() > 0 ? kChPos : kChNeg ;\r
122           Int_t id = fPID->GetParticleSpecie(partMC);\r
123           if(id != kSpUndefined) {\r
124             fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
125           }\r
126         }\r
127     }\r
128   \r
129   //main loop on tracks\r
130   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {\r
131     AliAODTrack* track = fAOD->GetTrack(iTracks);\r
132     if (!fTrackCuts->IsSelected(track)) continue;\r
133     \r
134     fPID->FillQAHistos(fHistMan, track, fTrackCuts);\r
135     \r
136     // //cut on q vectors track-by-track FIXME\r
137     // if ((qPos>fTrackCuts->GetQvecMin() && qPos<fTrackCuts->GetQvecMax() && track->Eta()<0) || (qNeg>fTrackCuts->GetQvecMin() && qNeg<fTrackCuts->GetQvecMax() && track->Eta()>=0)){\r
138     \r
139     //calculate DCA for AOD track\r
140     Double_t d[2], covd[3];\r
141     AliAODTrack* track_clone=(AliAODTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters\r
142     Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);\r
143     delete track_clone;\r
144     if(!isDCA)d[0]=-999;\r
145     \r
146     fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),d[0]);  // PT histo\r
147     \r
148     // get identity and charge\r
149     Int_t idRec  = fPID->GetParticleSpecie(track, fTrackCuts);\r
150     if(idRec == kSpUndefined) continue;\r
151     Int_t charge = track->Charge() > 0 ? kChPos : kChNeg;\r
152     \r
153     // rapidity cut\r
154     if(!fTrackCuts->CheckYCut ((AODParticleSpecies_t)idRec)) continue;\r
155     \r
156     // Fill histograms\r
157     fHistMan->GetHistogram2D(kHistPtRecSigma,idRec,charge)->Fill(track->Pt(),d[0]);\r
158     \r
159     /* MC Part */\r
160     if (arrayMC) {\r
161       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
162       if (!partMC) { \r
163         AliError("Cannot get MC particle");\r
164         continue; \r
165       }\r
166       // Check if it is primary, secondary from material or secondary from weak decay\r
167       Bool_t isPrimary           = partMC->IsPhysicalPrimary();\r
168       Bool_t isSecondaryMaterial = kFALSE; \r
169       Bool_t isSecondaryWeak     = kFALSE; \r
170       if(!isPrimary) {\r
171         Int_t mfl=-999,codemoth=-999;\r
172         Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()\r
173         if(indexMoth>=0){//is not fake\r
174           AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);\r
175           codemoth = TMath::Abs(moth->GetPdgCode());\r
176           mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));\r
177         }\r
178         if(mfl==3) isSecondaryWeak     = kTRUE; // add if(partMC->GetStatus() & kPDecay)? FIXME\r
179         else       isSecondaryMaterial = kTRUE;\r
180       }\r
181       \r
182      \r
183       // Get true ID\r
184       Int_t idGen     = fPID->GetParticleSpecie(partMC);\r
185       //if(TMath::Abs(partMC->Y())   > fTrackCuts->GetY()  ) continue;      // FIXME: do we need a rapidity cut on the generated?\r
186       // Fill histograms for primaries\r
187       if(idGen != kSpUndefined) {\r
188         \r
189         if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTrue,  idGen, charge)->Fill(track->Pt(),d[0]); \r
190         \r
191         if (isPrimary) {\r
192           fHistMan                    ->GetPtHistogram(kHistPtRecPrimary)->Fill(track->Pt(),d[0]);  // PT histo\r
193           fHistMan                    ->GetHistogram2D(kHistPtRecPrimary,      idGen, charge)->Fill(track->Pt(),d[0]);\r
194           if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTruePrimary,  idGen, charge)->Fill(track->Pt(),d[0]); \r
195           fHistMan                    ->GetHistogram2D(kHistPtRecSigmaPrimary, idRec, charge)->Fill(track->Pt(),d[0]); \r
196         }\r
197       }\r
198       //25th Apr - Muons are added to Pions -- FIXME\r
199       if ( partMC->PdgCode() == 13 && idRec == kSpPion) { \r
200         fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),d[0]); \r
201         if(isPrimary)\r
202           fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),d[0]); \r
203       }\r
204       if ( partMC->PdgCode() == -13 && idRec == kSpPion) { \r
205         fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),d[0]); \r
206         if (isPrimary) {\r
207           fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),d[0]); \r
208         }\r
209       }\r
210       \r
211       ///..... END FIXME\r
212       \r
213       // Fill secondaries\r
214       if(isSecondaryWeak    )  fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryWeakDecay, idRec, charge)->Fill(track->Pt(),d[0]);\r
215       if(isSecondaryMaterial)  fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryMaterial , idRec, charge)->Fill(track->Pt(),d[0]);\r
216       \r
217     }//end if(arrayMC)\r
218   } // end loop on tracks\r
219   \r
220   PostData(1, fHistMan  );\r
221   PostData(2, fEventCuts);\r
222   PostData(3, fTrackCuts);\r
223   PostData(4, fPID      );\r
224 }\r
225 \r
226 //_________________________________________________________________\r
227 void   AliAnalysisTaskSpectraAOD::Terminate(Option_t *)\r
228 {\r
229   // Terminate\r
230 }\r