TAmpt/macros/createAmptMC.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <Riostream.h>
   3 #include <TFile.h>
   4 #include <TTree.h>
   5 #include <TMath.h>
   6 #include <TParticle.h>
   7 #include <TPDGCode.h>
   8 #include <TDatabasePDG.h>
   9 #include <TRandom3.h>
  10 #include <TChain.h>
  11 #include <TFolder.h>
  12 #include "AliRun.h"
  13 #include "AliRunLoader.h"
  14 #include "AliHeader.h"
  15 #include "AliStack.h"
  16 #include "AliPDG.h"
  17 #include "AliGenAmpt.h"
  18 #include "TAmpt.h"
  19 #endif
  20
  21 class MyHeader
  22 {
  23 public:
  24   MyHeader() : fNATT(0), fEATT(0),fJATT(0),fNT(0),fNP(0),fN0(0),fN01(0),fN10(0),fN11(0),fBB(0),fRP(0) {;}
  25   virtual ~MyHeader() {;}
  26   Int_t       fNATT;
  27   Float_t     fEATT;
  28   Int_t       fJATT;
  29   Int_t       fNT;
  30   Int_t       fNP;
  31   Int_t       fN0;
  32   Int_t       fN01;
  33   Int_t       fN10;
  34   Int_t       fN11;
  35   Float_t     fBB;
  36   Float_t     fRP;
  37   ClassDef(MyHeader,1) // My header class
  38 };
  39
  40 class MyNuc : public TObject
  41 {
  42 public:
  43   MyNuc(Int_t pid=0, Int_t st=0, Int_t type=0, Double_t x=0, Double_t y=0, Double_t z=0) :
  44     TObject(), fPid(pid), fSt(st), fType(type), fX(x), fY(y), fZ(z) {;}
  45   Double_t    X()    const { return fX; }
  46   Double_t    Y()    const { return fY; }
  47   Double_t    Z()    const { return fZ; }
  48   Int_t       Pid()  const { return fPid; }
  49   Int_t       St()   const { return fSt; }
  50   Int_t       Type() const { return fType; }
  51 protected:
  52   Int_t       fPid;
  53   Int_t       fSt;
  54   Int_t       fType;
  55   Double32_t  fX;
  56   Double32_t  fY;
  57   Double32_t  fZ;
  58   ClassDef(MyNuc,1) // My nucleon class in cartesian coordinates
  59 };
  60
  61 class MyPart : public TObject
  62 {
  63 public:
  64   MyPart(Int_t pid=0, Int_t st=0, Int_t type=0, Double_t pt=0, Double_t eta=0, Double_t phi=0) :
  65     TObject(), fPid(pid), fSt(st), fType(type), fPt(pt), fEta(eta), fPhi(phi) {;}
  66   Double_t    Px()   const { return fPt*TMath::Cos(fPhi);  }
  67   Double_t    Py()   const { return fPt*TMath::Sin(fPhi);  }
  68   Double_t    Pz()   const { return fPt*TMath::SinH(fEta); }
  69   Double_t    Pt()   const { return fPt;  }
  70   Double_t    Eta()  const { return fEta; }
  71   Double_t    Phi()  const { return fPhi; }
  72   Int_t       Pid()  const { return fPid; }
  73   Int_t       St()   const { return fSt; }
  74   Int_t       Type() const { return fType; }
  75 protected:
  76   Int_t       fPid;
  77   Int_t       fSt;
  78   Int_t       fType;
  79   Double32_t  fPt;
  80   Double32_t  fEta;
  81   Double32_t  fPhi;
  82   ClassDef(MyPart,1) // My particle class in cylindrical coordinates
  83 };
  84
  85 void createAmptMC(Int_t nEvents,
  86                   const char *outFileName = "amptout.root");
  87
  88 void anaAmptMC(Int_t nEvents,
  89                const char *inFileNames = "/eliza6/alice/loizides/ampt/run_1/ampt*.root");
  90
  91 //-----------------------------------------------------------------------------------------------------
  92
  93 void createAmptMC(Int_t nEvents,
  94                   const char *outFileName)
  95 {
  96   TClass::GetClass("MyNuc")->IgnoreTObjectStreamer();
  97   TClass::GetClass("MyPart")->IgnoreTObjectStreamer();
  98
  99   AliPDG::AddParticlesToPdgDataBase();
 100   TDatabasePDG::Instance();
 101
 102   // Run loader
 103   TFolder *folder = new TFolder("myfolder","myfolder");
 104   AliRunLoader* rl = new AliRunLoader(folder);
 105   rl->MakeHeader();
 106   rl->MakeStack();
 107   AliStack* stack = rl->Stack();
 108   //AliHeader* rheader = rl->GetHeader();
 109
 110   AliGenAmpt *genHi = new AliGenAmpt(-1);
 111   genHi->SetStack(stack);
 112   genHi->SetEnergyCMS(2760);
 113   genHi->SetReferenceFrame("CMS");
 114   genHi->SetProjectile("A", 208, 82);
 115   genHi->SetTarget    ("A", 208, 82);
 116   genHi->SetPtHardMin (2);
 117   genHi->SetImpactParameterRange(0,30);
 118   genHi->SetJetQuenching(0); // enable jet quenching
 119   genHi->SetShadowing(1);    // enable shadowing
 120   genHi->SetDecaysOff(1);    // neutral pion and heavy particle decays switched off
 121   genHi->SetSpectators(0);   // track spectators
 122   genHi->SetIsoft(4);        // 4=string melting, 1=standard AMPT
 123   genHi->SetXmu(3.2264);     // parton xsection
 124   genHi->SetNtMax(150);      // time bins
 125   if (0) { //RHIC settings
 126     genHi->SetAlpha(0.47140452);
 127     genHi->SetStringFrag(2.2,0.5);
 128   }
 129   genHi->Init();
 130
 131   TClonesArray *inucs = new TClonesArray("TParticle",10000);
 132   TClonesArray *parts = new TClonesArray("MyPart",10000);
 133   TClonesArray *nucs  = new TClonesArray("MyNuc",500);
 134   MyHeader *myheader = new MyHeader;
 135
 136   TFile *outFile = TFile::Open(outFileName, "RECREATE");
 137   outFile->SetCompressionLevel(5);
 138   TDirectory::TContext context(outFile);
 139
 140   TTree *tree = new TTree("ampt", "AMPT tree");
 141   tree->Branch("parts",&parts, 256*1024, 99);
 142   tree->Branch("nucs",&nucs, 64*1024, 99);
 143   tree->Branch("info",&myheader, 32*1024, 99);
 144
 145   // create events and fill them
 146   for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {
 147
 148     cout << "Event " << iEvent+1 << "/" << nEvents << endl;;
 149     stack->Reset();
 150     genHi->Generate();
 151     parts->Clear();
 152     const TClonesArray *iarr = genHi->GetParticles();
 153     if (iarr) {
 154       for(Int_t i=0;i<iarr->GetEntriesFast();++i) {
 155         TParticle *p = (TParticle*)iarr->At(i);
 156         new((*parts)[i]) MyPart(p->GetPdgCode(),
 157                                 p->GetStatusCode(),
 158                                 p->GetUniqueID(),
 159                                 p->Pt(),
 160                                 p->Eta(),
 161                                 p->Phi());
 162       }
 163     }
 164     TAmpt *ampt = genHi->Ampt();
 165     if (ampt) {
 166       ampt->ImportNucleons(inucs);
 167       nucs->Clear();
 168       for(Int_t i=0;i<inucs->GetEntriesFast();++i) {
 169         TParticle *p = (TParticle*)inucs->At(i);
 170         new((*nucs)[i]) MyNuc(p->GetPdgCode(),
 171                               p->GetStatusCode(),
 172                               p->GetUniqueID(),
 173                               p->Vx(),
 174                               p->Vy(),
 175                               p->Vz());
 176       }
 177
 178       myheader->fNATT=ampt->GetNATT();
 179       myheader->fEATT=ampt->GetEATT();
 180       myheader->fJATT=ampt->GetJATT();
 181       myheader->fNT=ampt->GetNT();
 182       myheader->fNP=ampt->GetNP();
 183       myheader->fN0=ampt->GetN0();
 184       myheader->fN01=ampt->GetN01();
 185       myheader->fN10=ampt->GetN10();
 186       myheader->fN11=ampt->GetN11();
 187       myheader->fBB=ampt->GetBB();
 188       myheader->fRP=ampt->GetPhi();
 189     }
 190     tree->Fill();
 191   } // end of event loop
 192   tree->Write();
 193   outFile->Close();
 194 }
 195
 196 //-----------------------------------------------------------------------------------------------------
 197
 198 void anaAmptMC(Int_t nEvents,
 199                const char *inFileNames)
 200 {
 201
 202   TChain *c = new TChain("ampt");
 203   c->Add(inFileNames);
 204
 205
 206   TClonesArray *parts = 0;
 207   TClonesArray *nucs  = 0;
 208   MyHeader *info = 0;
 209
 210   c->SetBranchAddress("info",&info);
 211   c->SetBranchAddress("nucs",&nucs);
 212   c->SetBranchAddress("parts",&parts);
 213
 214   Int_t nRead = nEvents;
 215   if (nRead<0)
 216     nRead = c->GetEntries();
 217   else if (0 && (nRead>c->GetEntries()))
 218     nRead = c->GetEntries();
 219
 220    for (Int_t ev=0;ev<nRead;++ev) {
 221       c->GetEntry(ev);
 222
 223       Int_t fAN=0, fBN=0, fAStat[250], fBStat[250];
 224       Double_t fXA[250], fXB[250], fYA[250], fYB[250];
 225       for (Int_t k=0;k<nucs->GetEntries();++k) {
 226         MyNuc *n = (MyNuc*)nucs->At(k);
 227         if (n->Type()>0) {
 228           fAStat[fAN] =0;
 229           fXA[fAN] = n->X();
 230           fYA[fAN] = n->Y();
 231           ++fAN;
 232         } else {
 233           fBStat[fBN] =0;
 234           fXB[fBN] = n->X();
 235           fYB[fBN] = n->Y();
 236           ++fBN;
 237         }
 238       }
 239
 240       // Glauber calculation
 241       Double_t fXSect = 65; //mbarn
 242       Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
 243
 244       // For each of the A nucleons in nucleus B
 245       for (Int_t i = 0; i<fBN; i++) {
 246         for (Int_t j = 0 ; j < fAN ;j++) {
 247           Double_t dx = fXB[i]-fXA[j];
 248           Double_t dy = fYB[i]-fYA[j];
 249           Double_t dij = dx*dx+dy*dy;
 250           if (dij < d2) {
 251             fBStat[i]++;
 252             fAStat[j]++;
 253           }
 254         }
 255       }
 256       // Calculate npart
 257       Int_t npart=0;
 258       for (Int_t i = 0; i<fAN; i++) {
 259         if (fAStat[i]>0) {
 260           npart++;
 261         }
 262       }
 263       for (Int_t i = 0; i<fBN; i++) {
 264         if (fBStat[i]>0) {
 265           npart++;
 266         }
 267       }
 268       cout << ev << " : "  << info->fBB << " np " << npart << " vs " << info->fNP+info->fNT << endl;
 269    }
 270 }