MUON/MUONdigits.C

   1 void MUONdigits (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nCathode=1)
   2 {
   3 /////////////////////////////////////////////////////////////////////////
   4 //   This macro is a small example of a ROOT macro
   5 //   illustrating how to read the output of GALICE
   6 //   and do some analysis.
   7 //
   8 /////////////////////////////////////////////////////////////////////////
   9
  10 // Dynamically link some shared libs
  11
  12    if (gClassTable->GetID("AliRun") < 0) {
  13        gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so");
  14        gSystem->Load("$ALITOP/cern.so/lib/libPythia.so");
  15        gSystem->Load("$ROOTSYS/lib/libEG.so");
  16        gSystem->Load("$ROOTSYS/lib/libEGPythia.so");
  17        gSystem->Load("libGeant3Dummy.so");        //a dummy version of Geant3
  18        gSystem->Load("PHOS/libPHOSdummy.so");     //the standard Alice classes
  19        gSystem->Load("libgalice.so");             //the standard Alice classes
  20    }
  21
  22
  23 // Connect the Root Galice file containing Geometry, Kine and Hits
  24
  25    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  26    if (file) file->Close();
  27    file = new TFile("galice.root","UPDATE");
  28    file->ls();
  29    //   file->Map();
  30
  31    printf ("I'm after Map \n");
  32
  33 // Get AliRun object from file or create it if not on file
  34
  35    if (!gAlice) {
  36        gAlice = (AliRun*)file->Get("gAlice");
  37        if (gAlice) printf("AliRun object found on file\n");
  38        if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  39    }
  40    printf ("I'm after gAlice \n");
  41
  42    AliMUON *MUON  = gAlice->GetDetector("MUON");
  43
  44    AliMUONchamber*  iChamber;
  45    AliMUONsegmentation*  segmentation;
  46
  47    Int_t Npx[10];
  48    Int_t Npy[10];
  49    Int_t trk[50];
  50    Int_t chtrk[50];
  51
  52
  53    Int_t nxmax=1026;
  54    Int_t nymax=1026;
  55    AliMUONlist *elem[10*1026*1026];
  56    AliMUONlist **ppe=elem;
  57    TObjArray *obj=new TObjArray;
  58
  59    Int_t digits[3];
  60    // fill the info array
  61    TVector *trinfo;
  62    trinfo=new TVector(2);
  63
  64 //
  65 //   Loop over events
  66 //
  67
  68    Int_t Nh=0;
  69    Int_t Nh1=0;
  70    for (int nev=0; nev<= evNumber2; nev++) {
  71        Int_t nparticles = gAlice->GetEvent(nev);
  72        cout << "nev         " <<nev<<endl;
  73        cout << "nparticles  " <<nparticles<<endl;
  74        if (nev < evNumber1) continue;
  75        if (nparticles <= 0) return;
  76
  77        TTree *TH = gAlice->TreeH();
  78        Int_t ntracks = TH->GetEntries();
  79        Int_t Nc=0;
  80        //
  81        Int_t counter=0;
  82
  83        // loop over cathodes
  84        for (int icat=0;icat<nCathode;icat++) {
  85
  86            //initialize the arrray of pointers !!!!!
  87            ppe=memset(elem,0,sizeof(void*)*10*1026*1026);
  88
  89            printf("Start loop over tracks \n");
  90 //
  91 //   Loop over events
  92 //
  93            for (Int_t track=0; track<ntracks;track++) {
  94                gAlice->ResetHits();
  95                Int_t nbytes += TH->GetEvent(track);
  96                if (MUON)  {
  97                    for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
  98                        mHit;
  99                        mHit=(AliMUONhit*)MUON->NextHit())
 100                    {
 101                        Int_t   nch   = mHit->fChamber;  // chamber number
 102                        Float_t x     = mHit->fX;        // x-pos of hit
 103                        Float_t y     = mHit->fY;        // y-pos
 104 //
 105 //
 106                        if (nch >10) continue;
 107
 108                        iChamber = &(MUON->Chamber(nch-1));
 109                        response=iChamber->GetResponseModel();
 110
 111                        Int_t nsec=iChamber->Nsec();
 112                        Int_t rmin = (Int_t)iChamber->frMin;
 113                        Int_t rmax = (Int_t)iChamber->frMax;
 114 //
 115 //
 116                        for (AliMUONcluster* mPad=(AliMUONcluster*)MUON->FirstPad(mHit);
 117                             mPad;
 118                             mPad=(AliMUONcluster*)MUON->NextPad())
 119                        {
 120                            Int_t cathode     = mPad->fCathode;   // chamber number
 121                            Int_t nhit     = mPad->fHitNumber; // hit number
 122                            Int_t qtot     = mPad->fQ;         // charge
 123                            Int_t ipx      = mPad->fPadX;      // pad number on X
 124                            Int_t ipy      = mPad->fPadY;      // pad number on Y
 125                            Int_t iqpad    = mPad->fQpad;      // charge per pad
 126                            Int_t izone    = mPad->fRSec;      // r-pos of pad
 127 //
 128 //
 129                            if (cathode != (icat+1)) continue;
 130                            segmentation=iChamber->GetSegmentationModel(cathode);
 131                            Int_t Npx[nch-1]  = segmentation->Npx();
 132                            Int_t Npy[nch-1]  = segmentation->Npy();
 133
 134                            Int_t npx=Npx[nch-1];
 135                            Int_t npy=Npy[nch-1];
 136                            Float_t thex, they;
 137
 138                            segmentation->GetPadCxy(ipx,ipy,thex,they);
 139                            Float_t rpad=TMath::Sqrt(thex*thex+they*they);
 140
 141                            if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
 142
 143
 144                            trinfo(0)=(Float_t)track;
 145                            trinfo(1)=(Float_t)iqpad;
 146
 147                            digits[0]=ipx;
 148                            digits[1]=ipy;
 149                            digits[2]=iqpad;
 150
 151                            // build the list of fired pads and update the info
 152                            AliMUONlist *pdigit=
 153                                elem[(nch-1)*(1026*1026)+(ipy+npy)*1026+(ipx+npx)];
 154                            if (pdigit==0) {
 155                                obj->AddAtAndExpand(new AliMUONlist(nch-1,digits),counter);
 156                                counter++;
 157                                Int_t last=obj->GetLast();
 158                                pdigit=(AliMUONlist*)obj->At(last);
 159                                elem[(nch-1)*(1026*1026)+(ipy+npy)*1026+(ipx+npx)]=pdigit;
 160                                // list of tracks
 161                                TObjArray *trlist=(TObjArray*)pdigit->TrackList();
 162                                trlist->Add(trinfo);
 163                            } else {
 164                                // update charge
 165                                (*pdigit).fSignal+=iqpad;
 166                                // update list of tracks
 167                                TObjArray* trlist=(TObjArray*)pdigit->TrackList();
 168                                Int_t last_entry=trlist->GetLast();
 169                                TVector *ptrk=(TVector*)trlist->At(last_entry);
 170                                Int_t last_track=Int_t(ptrk(0));
 171                                Int_t last_charge=Int_t(ptrk(1));
 172                                if (last_track==track) {
 173                                    last_charge+=iqpad;
 174                                    trlist->RemoveAt(last_entry);
 175                                    trinfo(0)=last_track;
 176                                    trinfo(1)=last_charge;
 177                                    trlist->AddAt(trinfo,last_entry);
 178                                } else {
 179                                    trlist->Add(trinfo);
 180                                }
 181                                // check the track list
 182                                Int_t nptracks=trlist->GetEntriesFast();
 183                                if (nptracks > 2) {
 184                                    printf("Attention - nptracks > 2  %d \n",nptracks);
 185                                    printf("cat,nch,ix,iy %d %d %d %d  \n",icat+1,nch,ipx,ipy);
 186                                    for (Int_t tr=0;tr<nptracks;tr++) {
 187                                        TVector *pptrk=(TVector*)trlist->At(tr);
 188                                        trk[tr]=Int_t(pptrk(0));
 189                                        chtrk[tr]=Int_t(pptrk(1));
 190                                    }
 191                                } // end if nptracks
 192                            } //  end if pdigit
 193                        } //end loop over clusters
 194                    } // hit loop
 195                } // if MUON
 196            } // track loop
 197
 198            Int_t tracks[10];
 199            Int_t charges[10];
 200            cout<<"start filling digits "<<endl;
 201            const Int_t zero_supm = 6;
 202            const Int_t adc_satm = 1024;
 203            Int_t nd=0;
 204
 205
 206            Int_t nentries=obj->GetEntriesFast();
 207            printf(" nentries %d \n",nentries);
 208
 209            // start filling the digits
 210
 211            for (Int_t nent=0;nent<nentries;nent++) {
 212                AliMUONlist *address=(AliMUONlist*)obj->At(nent);
 213                if (address==0) continue;
 214                // do zero-suppression and signal truncation
 215                Int_t ich=address->fRpad;
 216                Int_t q=address->fSignal;
 217                // if ( q <= zero_supm  || rpad > 55) continue;
 218                if ( q <= zero_supm ) continue;
 219                if ( q > adc_satm)  q=adc_satm;
 220                digits[0]=address->fPadX;
 221                digits[1]=address->fPadY;
 222                digits[2]=q;
 223
 224                TObjArray* trlist=(TObjArray*)address->TrackList();
 225                Int_t nptracks=trlist->GetEntriesFast();
 226                // this was changed to accomodate the real number of tracks
 227                if (nptracks > 10) {
 228                    cout<<"Attention - nptracks > 3 "<<nptracks<<endl;
 229                    nptracks=10;
 230                }
 231                if (nptracks > 2) {
 232                    printf("Attention - nptracks > 2  %d \n",nptracks);
 233                    printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
 234                }
 235                for (Int_t tr=0;tr<nptracks;tr++) {
 236                    TVector *pp=(TVector*)trlist->At(tr);
 237                    Int_t entries=pp->GetNrows();
 238                    tracks[tr]=Int_t(pp(0));
 239                    charges[tr]=Int_t(pp(1));
 240                }      //end loop over list of tracks for one pad
 241
 242                // fill digits
 243                MUON->AddDigits(ich,tracks,charges,digits);
 244            }
 245            cout<<"I'm out of the loops for digitisation"<<endl;
 246            gAlice->TreeD()->Fill();
 247            TTree *TD=gAlice->TreeD();
 248            int ndig=TD->GetEntries();
 249            cout<<"number of digits  "<<ndig<<endl;
 250            TClonesArray *fDch;
 251            for (int i=0;i<10;i++) {
 252                fDch= MUON->DigitsAddress(i);
 253                int ndig=fDch->GetEntriesFast();
 254                printf (" i, ndig %d %d \n",i,ndig);
 255            }
 256            MUON->ResetDigits();
 257            obj->Clear();
 258        } //end loop over cathodes
 259        char hname[30];
 260        sprintf(hname,"TreeD%d",nev);
 261        gAlice->TreeD()->Write(hname);
 262        file->ls();
 263    } // event loop
 264    file->Close();
 265 //   cout<<"END  digitisation "<<endl;
 266 }