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