/************************************************************************** * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ //----------------------------------------------------------------- // AliAnalysisTaskSpectraAOD class //----------------------------------------------------------------- #include "TChain.h" #include "TTree.h" #include "TLegend.h" #include "TH1F.h" #include "TH2F.h" #include "TCanvas.h" #include "AliAnalysisTask.h" #include "AliAnalysisManager.h" #include "AliAODTrack.h" #include "AliAODMCParticle.h" #include "AliVParticle.h" #include "AliAODEvent.h" #include "AliAODInputHandler.h" #include "AliAnalysisTaskSpectraAOD.h" #include "AliAnalysisTaskESDfilter.h" #include "AliAnalysisDataContainer.h" #include "AliSpectraAODHistoManager.h" #include "AliSpectraAODTrackCuts.h" #include "AliSpectraAODEventCuts.h" #include "AliCentrality.h" #include "TProof.h" #include "AliPID.h" #include "AliVEvent.h" #include "AliPIDResponse.h" #include "AliStack.h" #include "AliSpectraAODPID.h" #include #include using namespace AliSpectraNameSpace; using namespace std; ClassImp(AliAnalysisTaskSpectraAOD) // EX1 This stuff tells root to implement the streamer, inspector methods etc (we discussed about it today) //________________________________________________________________________ AliAnalysisTaskSpectraAOD::AliAnalysisTaskSpectraAOD(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0), fPID(0), fIsMC(0), fNRebin(0) { // Default constructor DefineInput(0, TChain::Class()); DefineOutput(1, AliSpectraAODHistoManager::Class()); DefineOutput(2, AliSpectraAODEventCuts::Class()); DefineOutput(3, AliSpectraAODTrackCuts::Class()); DefineOutput(4, AliSpectraAODPID::Class()); fNRebin=0; } //________________________________________________________________________ //________________________________________________________________________ void AliAnalysisTaskSpectraAOD::UserCreateOutputObjects() { // create output objects fHistMan = new AliSpectraAODHistoManager("SpectraHistos",fNRebin); if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro"); if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro"); if (!fPID) AliFatal("PID object should be set in the steering macro"); PostData(1, fHistMan ); PostData(2, fEventCuts); PostData(3, fTrackCuts); PostData(4, fPID ); } //________________________________________________________________________ void AliAnalysisTaskSpectraAOD::UserExec(Option_t *) { // main event loop //Printf("ALIVE"); fAOD = dynamic_cast(fInputEvent); if (!fAOD) { AliWarning("ERROR: AliAODEvent not available \n"); return; } if (strcmp(fAOD->ClassName(), "AliAODEvent")) { AliFatal("Not processing AODs"); } if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection //AliCentrality fAliCentral*; // if ((fAOD->GetCentrality())->GetCentralityPercentile("V0M") > 5.) return; // First do MC to fill up the MC particle array, such that we can use it later TClonesArray *arrayMC = 0; if (fIsMC) { arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName()); if (!arrayMC) { AliFatal("Error: MC particles branch not found!\n"); } Int_t nMC = arrayMC->GetEntries(); for (Int_t iMC = 0; iMC < nMC; iMC++) { AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC); if(!partMC->Charge()) continue;//Skip neutrals //if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax()){//charged hadron are filled inside the eta acceptance //Printf("%f %f-%f",partMC->Eta(),fTrackCuts->GetEtaMin(),fTrackCuts->GetEtaMax()); if(partMC->Eta() < fTrackCuts->GetEtaMin() || partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!! fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); //rapidity cut if(TMath::Abs(partMC->Y()) > fTrackCuts->GetY() ) continue; // check for true PID + and fill P_t histos Int_t charge = partMC->Charge() > 0 ? kChPos : kChNeg ; Int_t id = fPID->GetParticleSpecie(partMC); if(id != kSpUndefined) { fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary()); } } } Double_t mass[3]={1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};//FIXME masses to be taken from AliHelperPID //main loop on tracks for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) { AliAODTrack* track = dynamic_cast(fAOD->GetTrack(iTracks)); if(!track) AliFatal("Not a standard AOD"); if (!fTrackCuts->IsSelected(track,kTRUE)) continue; fPID->FillQAHistos(fHistMan, track, fTrackCuts); //calculate DCA for AOD track Double_t dca=track->DCA(); if(dca==-999.){// track->DCA() does not work in old AOD production Double_t d[2], covd[3]; AliAODTrack* track_clone=(AliAODTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd); delete track_clone; if(!isDCA)d[0]=-999.; dca=d[0]; } fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),dca); // PT histo // get identity and charge Int_t idRec = fPID->GetParticleSpecie(fHistMan,track, fTrackCuts); Int_t charge = track->Charge() > 0 ? kChPos : kChNeg; // Fill histograms, only if inside y and nsigma acceptance if(idRec != kSpUndefined){ if(fTrackCuts->CheckYCut (mass[idRec]))fHistMan->GetHistogram2D(kHistPtRecSigma,idRec,charge)->Fill(track->Pt(),dca); }//can't put a continue because we still have to fill allcharged primaries, done later /* MC Part */ if (arrayMC) { AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel())); if (!partMC) { AliError("Cannot get MC particle"); continue; } // Check if it is primary, secondary from material or secondary from weak decay Bool_t isPrimary = partMC->IsPhysicalPrimary(); Bool_t isSecondaryMaterial = kFALSE; Bool_t isSecondaryWeak = kFALSE; if(!isPrimary) { Int_t mfl=-999,codemoth=-999; Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother() if(indexMoth>=0){//is not fake AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth); codemoth = TMath::Abs(moth->GetPdgCode()); mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth)))); } //Int_t uniqueID = partMC->GetUniqueID(); //cout<<"uniqueID: "<GetUniqueID()<<" "<GetStatus()<<" "<GetStatus() & kPDecay)? FIXME else isSecondaryMaterial = kTRUE; } if (isPrimary)fHistMan->GetPtHistogram(kHistPtRecPrimary)->Fill(track->Pt(),dca); // PT histo of primaries //nsigma cut (reconstructed nsigma) if(idRec == kSpUndefined) continue; // rapidity cut (reconstructed pt and identity) if(!fTrackCuts->CheckYCut (mass[idRec])) continue; // Get true ID Int_t idGen = fPID->GetParticleSpecie(partMC); //if(TMath::Abs(partMC->Y()) > fTrackCuts->GetY() ) continue; // FIXME: do we need a rapidity cut on the generated? // Fill histograms for primaries if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTrue, idGen, charge)->Fill(track->Pt(),dca); if (isPrimary) { fHistMan ->GetHistogram2D(kHistPtRecSigmaPrimary, idRec, charge)->Fill(track->Pt(),dca); if(idGen != kSpUndefined) { fHistMan ->GetHistogram2D(kHistPtRecPrimary, idGen, charge)->Fill(track->Pt(),dca); if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTruePrimary, idGen, charge)->Fill(track->Pt(),dca); } } //25th Apr - Muons are added to Pions -- FIXME if ( partMC->PdgCode() == 13 && idRec == kSpPion) { fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),dca); if(isPrimary) fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),dca); } if ( partMC->PdgCode() == -13 && idRec == kSpPion) { fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),dca); if (isPrimary) { fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),dca); } } ///..... END FIXME // Fill secondaries if(isSecondaryWeak ) fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryWeakDecay, idRec, charge)->Fill(track->Pt(),dca); if(isSecondaryMaterial) fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryMaterial , idRec, charge)->Fill(track->Pt(),dca); }//end if(arrayMC) } // end loop on tracks PostData(1, fHistMan ); PostData(2, fEventCuts); PostData(3, fTrackCuts); PostData(4, fPID ); } //_________________________________________________________________ void AliAnalysisTaskSpectraAOD::Terminate(Option_t *) { // Terminate }