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))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;\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       Bool_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 \r
 187       // Fill histograms for primaries\r
 188       if (isPrimary) {\r
 189         fHistMan                    ->GetHistogram2D(kHistPtRecPrimary,      idGen, charge)->Fill(track->Pt(),d[0]);  // PT histo\r
 190         if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTruePrimary,  idGen, charge)->Fill(track->Pt(),d[0]); \r
 191         fHistMan                    ->GetHistogram2D(kHistPtRecSigmaPrimary, idRec, charge)->Fill(track->Pt(),d[0]); \r
 192               \r
 193       }\r
 194       //25th Apr - Muons are added to Pions -- FIXME\r
 195       if ( partMC->PdgCode() == 13 && idRec == kSpPion) { \r
 196         fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),d[0]); \r
 197         if(isPrimary)\r
 198           fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),d[0]); \r
 199       }\r
 200       if ( partMC->PdgCode() == -13 && idRec == kSpPion) { \r
 201         fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),d[0]); \r
 202         if (isPrimary) {\r
 203           fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),d[0]); \r
 204         }\r
 205       }\r
 206 \r
 207       ///..... END FIXME\r
 208 \r
 209       // Fill secondaries\r
 210       if(isSecondaryWeak    )  fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryWeakDecay, idRec, charge)->Fill(track->Pt(),d[0]);\r
 211       if(isSecondaryMaterial)  fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryMaterial , idRec, charge)->Fill(track->Pt(),d[0]);\r
 212             \r
 213     }//end if(arrayMC)\r
 214   } // end loop on tracks\r
 215   \r
 216   PostData(1, fHistMan);\r
 217   PostData(2, fEventCuts);\r
 218   PostData(3, fTrackCuts);\r
 219   PostData(4, fTrackCuts);\r
 220 }\r
 221 \r
 222 //_________________________________________________________________\r
 223 void   AliAnalysisTaskSpectraAOD::Terminate(Option_t *)\r
 224 {\r
 225   // Terminate\r
 226 }\r