MUON/AliMUONproto.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /*
  19 Implementation of MUON Chamber Prototype Class
  20 Both read digits from raw data or use the Monte-Carlo.
  21 1-Feb-2000 Rachid GUERNANE, IPN Lyon guernane@ipnl.in2p3.fr
  22 */
  23
  24
  25 ////////////////////////////////////////////////
  26 //  Manager and hits classes for set:PROTO    //
  27 ////////////////////////////////////////////////
  28
  29 #include <TTUBE.h>
  30 #include <TBRIK.h>
  31 #include <TRotMatrix.h>
  32 #include <TNode.h>
  33 #include <TTree.h>
  34 #include <TRandom.h>
  35 #include <TObject.h>
  36 #include <TVector.h>
  37 #include <TObjArray.h>
  38 #include <TMinuit.h>
  39 #include <TParticle.h>
  40 #include <TROOT.h>
  41 #include <TFile.h>
  42 #include <TNtuple.h>
  43 #include <TCanvas.h>
  44 #include <TPad.h>
  45 #include <TDirectory.h>
  46 #include <TObjectTable.h>
  47 #include <AliPDG.h>
  48 #include <TArrayI.h>
  49 #include <AliDetector.h>
  50
  51 #include "AliMUONChamber.h"
  52 #include "AliMUONproto.h"
  53 #include "AliMUONHit.h"
  54 #include "TTUBE.h"
  55 #include "AliMUONClusterFinder.h"
  56 #include "AliRun.h"
  57 #include "Riostream.h"
  58 #include "AliCallf77.h"
  59 #include "AliConst.h"
  60 //#include "chainalice2.h"
  61 #include "AliMUONSegmentationV0.h"
  62 //#include "AliMUONSegResV11.h"
  63 #include "AliLog.h"
  64
  65 ClassImp(AliMUONproto)
  66
  67 //___________________________________________
  68 AliMUONproto::AliMUONproto()
  69         : AliMUON()
  70 {
  71     AliInfo(" Calling AliMUONproto constructor...");
  72
  73         //
  74         //
  75         //
  76 }
  77
  78 //___________________________________________
  79 AliMUONproto::AliMUONproto(const char *name, const char *title)
  80         : AliMUON(name,title)
  81 //  AliMUON defines everything, but the chamber for NCH=1
  82 {
  83
  84 //  z-Positions of Chambers
  85     const Float_t zch[1] = {975.};
  86 //
  87 //  inner diameter
  88     const Float_t dmi[1] = {0.};
  89 //
  90 //  outer diameter
  91     const Float_t dma[1] = {29.2};
  92 //
  93 //
  94 //    Default Parameters for ALICE2 prototype
  95
  96 //
  97     (*fChambers)[0] = new AliMUONChamber();
  98     AliMUONChamber* chamber = (AliMUONChamber*) (*fChambers)[0];
  99     chamber->SetGid(0);
 100     chamber->SetZ(zch[0]);
 101 //
 102     chamber->InitGeo(zch[0]);
 103     chamber->SetRInner(dmi[0]/2);
 104     chamber->SetROuter(dma[0]/2);
 105
 106     for (int i = 0; i <= 99; i++) {
 107         fThreshold[i] = 0.;
 108     }
 109
 110 }
 111
 112 #ifdef WE_FORGET_FOR_THE_MOMENT
 113 //___________________________________________
 114 void AliMUONproto::GetRawDigits(Int_t evnb, Int_t *lptr, Int_t ilen) {
 115
 116     Int_t ip = 0;
 117     Int_t equip = 0;
 118     Int_t nochnl;
 119     Int_t loop;
 120     Int_t val;
 121     Int_t itype;
 122     Int_t id;
 123     Int_t serial;
 124     Int_t equiplen;
 125
 126     Int_t digits[5];
 127
 128     AliMUON* MUON  = (AliMUON*)gAlice->GetModule("MUON");
 129     AliMUONSegmentationV0* seg = (AliMUONSegmentationV0*) Chamber(0).SegmentationModel(1);
 130
 131     Int_t tracks[10];
 132     Int_t charges[10];
 133
 134     for (Int_t i = 0; i < 10; i++) {
 135         tracks[i] = 0;
 136         charges[i] = 0;
 137     }
 138
 139     Int_t ich = 0;
 140
 141   nwtype:
 142
 143     itype = lptr[ip++];
 144     id = lptr[ip++];
 145     equiplen = lptr[ip++];
 146
 147     if (equiplen < 0 ) return;
 148
 149     if (itype != (int)(0XCACCA008)) {
 150         ip += equiplen;
 151         if (ip < ilen-2)  goto nwtype;
 152     }
 153     else {
 154         serial = id >> 16;
 155         equip = id & 0x1;
 156         if ((serial == 190) && (equip == 1)) {
 157             for (loop = 0; loop < equiplen; loop++) {
 158                 nochnl = (lptr[ip] & 0x7ff000 ) >> 12;
 159                 val = lptr[ip] & 0x3ff;
 160                 // fill digits from raw data according to cathode connexions
 161                 if (group[nochnl][2][1]!=0)
 162                     digits[0] = group[nochnl][2][1] - 12;
 163                 else if  (group[nochnl][1][1]!=0)
 164                     digits[0] = group[nochnl][1][1] - 12;
 165                 else
 166                     digits[0] = group[nochnl][0][1] - 12;
 167                 digits[1] = group[nochnl][0][0];
 168                 if (digits[0] != seg->Ix(digits[0], digits[1]))
 169                     printf("Pb pour ix=%d,iy=%d\n", digits[0], digits[1]);
 170                 digits[2] = val;
 171                 digits[3] = 0;
 172                 digits[4] = 0;
 173                 if (digits[2] >= fThreshold[nochnl])
 174                     MUON->AddDigits(ich, tracks, charges, digits);
 175
 176                 ip++;
 177             }
 178         }
 179         else
 180             ip += equiplen;
 181
 182         if (ip < ilen-2) goto nwtype;
 183     }
 184
 185     gAlice->TreeD()->Fill();
 186     MUON->ResetDigits();
 187
 188     gAlice->TreeD()->Fill();
 189     MUON->ResetDigits();
 190
 191     char hname[30];
 192     sprintf(hname, "TreeD%d", evnb);
 193     gAlice->TreeD()->Write(hname,TObject::kOverwrite);
 194     // reset tree
 195     gAlice->TreeD()->Reset();
 196
 197 }
 198 #endif
 199
 200 //___________________________________________
 201 void AliMUONproto::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
 202 {
 203     Int_t i=2*(id-1);
 204     ((AliMUONChamber*)(*fChambers)[i])->SetPadSize(isec,p1,p2);
 205 }
 206
 207 //___________________________________________
 208 void AliMUONproto::SetChargeSlope(Int_t id, Float_t p1)
 209 {
 210     Int_t i=2*(id-1);
 211     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSlope(p1);
 212 }
 213
 214 //___________________________________________
 215 void AliMUONproto::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
 216 {
 217     Int_t i=2*(id-1);
 218     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
 219 }
 220
 221 //___________________________________________
 222 void AliMUONproto::SetMaxAdc(Int_t id, Float_t p1)
 223 {
 224     Int_t i=2*(id-1);
 225     ((AliMUONChamber*) (*fChambers)[i])->SetMaxAdc(p1);
 226 }
 227
 228 //___________________________________________
 229 void AliMUONproto::CreateGeometry()
 230 {
 231     Int_t *idtmed = fIdtmed->GetArray()-1099;
 232 //
 233     Float_t zpos;
 234     Float_t bpar[3];
 235     Int_t idGas=idtmed[1100];
 236
 237     AliMUONChamber *iChamber;
 238 //********************************************************************
 239 //                          Prototype ALICE2                        **
 240 //********************************************************************
 241     iChamber=(AliMUONChamber*) (*fChambers)[0];
 242
 243     zpos=iChamber->Z();
 244 //
 245     const Float_t X_POS = 11*.975/2; //half x
 246     const Float_t Y_POS = 18*.55/2;  //half y
 247     const Float_t Z_POS = 0.325;
 248
 249     bpar[0] = X_POS;
 250     bpar[1] = Y_POS;
 251     bpar[2] = Z_POS;
 252
 253     gMC->Gsvolu("C01G", "BOX", idGas, bpar, 3);
 254     gMC->Gspos("C01G", 1, "ALIC", -bpar[0], bpar[1], zpos, 0, "ONLY");
 255 }
 256
 257 //___________________________________________
 258 void AliMUONproto::CreateMaterials()
 259 {
 260   // *** DEFINITION OF AVAILABLE MUON MATERIALS ***
 261   //
 262   //     Ar-CO2 gas
 263     Float_t ag1[3]   = { 39.95,12.01,16. };
 264     Float_t zg1[3]   = { 18.,6.,8. };
 265     Float_t wg1[3]   = { .8,.0667,.13333 };
 266     Float_t dg1      = .001821;
 267
 268     Float_t epsil  = .001; // Tracking precision,
 269     Float_t tmaxfd = -20.; // Maximum angle due to field deflection
 270     Float_t stmin  = -.8;
 271     Int_t ISXFLD   = gAlice->Field()->Integ();
 272     Float_t SXMGMX = gAlice->Field()->Max();
 273
 274     //
 275     // --- Define the various materials for GEANT ---
 276     AliMixture(22, "ArCO2 80%$", ag1, zg1, dg1, 3, wg1);
 277
 278     //
 279     //    Ar/CO2
 280     AliMedium(1, "ARG_CO2   ", 22, 1, ISXFLD, SXMGMX, tmaxfd, fMaxStepGas,
 281               fMaxDestepAlu, epsil, stmin);
 282     //    Air
 283         //AliMedium(1, "AIR_CH_US         ", 15, 1, ISXFLD, SXMGMX, tmaxfd, -1., -.3, epsil, stmin);
 284 }
 285
 286 //___________________________________________
 287 void AliMUONproto::Init()
 288 {
 289    printf("\n\n\n Start Init for Prototype ALICE2 - CPC chamber type\n\n\n");
 290
 291    //
 292    // Initialize Tracking Chambers
 293    //
 294
 295    ( (AliMUONChamber*) (*fChambers)[0])->Init();
 296
 297    //
 298    // Set the chamber (sensitive region) GEANT identifier
 299    ((AliMUONChamber*)(*fChambers)[0])->SetGid(gMC->VolId("C01G"));
 300
 301    printf("\n\n\n Finished Init for Prototype ALICE2 - CPC chamber type\n\n\n");
 302
 303 }
 304
 305 //___________________________________________
 306 void AliMUONproto::StepManager()
 307 {
 308   Int_t          copy, id;
 309   static Int_t   idvol;
 310   static Int_t   vol[2];
 311   Int_t          ipart;
 312   TLorentzVector pos;
 313   TLorentzVector mom;
 314   Float_t        theta,phi;
 315   Float_t        destep, step;
 316
 317   static Float_t eloss, eloss2, xhit, yhit, tlength;
 318   const  Float_t big=1.e10;
 319
 320   //  modifs perso
 321   static Float_t hits[14];
 322
 323   TClonesArray &lhits = *fHits;
 324
 325   //
 326   // Set maximum step size for gas
 327   // numed=gMC->GetMedium();
 328   //
 329   // Only charged tracks
 330   if( !(gMC->TrackCharge()) ) return;
 331   //
 332   // Only gas gap inside chamber
 333   // Tag chambers and record hits when track enters
 334   idvol=-1;
 335   id=gMC->CurrentVolID(copy);
 336
 337     for (Int_t i=1; i<=kNCH; i++) {
 338       if(id==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()){
 339           vol[0]=i;
 340           idvol=i-1;
 341       }
 342     }
 343     if (idvol == -1) return;
 344   //
 345   // Get current particle id (ipart), track position (pos)  and momentum (mom)
 346   gMC->TrackPosition(pos);
 347   gMC->TrackMomentum(mom);
 348
 349   ipart  = gMC->TrackPid();
 350   //Int_t ipart1 = gMC->IdFromPDG(ipart);
 351   //printf("ich, ipart %d %d \n",vol[0],ipart1);
 352
 353   //
 354   // momentum loss and steplength in last step
 355   destep = gMC->Edep();
 356   step   = gMC->TrackStep();
 357
 358   //
 359   // record hits when track enters ...
 360   if( gMC->IsTrackEntering()) {
 361       gMC->SetMaxStep(fMaxStepGas);
 362       Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 363       Double_t rt = TMath::Sqrt(tc);
 364       Double_t pmom = TMath::Sqrt(tc+mom[2]*mom[2]);
 365       Double_t tx=mom[0]/pmom;
 366       Double_t ty=mom[1]/pmom;
 367       Double_t tz=mom[2]/pmom;
 368       Double_t s=((AliMUONChamber*)(*fChambers)[idvol])
 369           ->ResponseModel()
 370           ->Pitch()/tz;
 371       theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 372       phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 373       hits[0] = Float_t(ipart);         // Geant3 particle type
 374       hits[1] = pos[0]+s*tx;                 // X-position for hit
 375       hits[2] = pos[1]+s*ty;                 // Y-position for hit
 376       hits[3] = pos[2]+s*tz;                 // Z-position for hit
 377       hits[4] = theta;                  // theta angle of incidence
 378       hits[5] = phi;                    // phi angle of incidence
 379       hits[8] = (Float_t) fNPadHits;   // first padhit
 380       hits[9] = -1;                     // last pad hit
 381
 382       // modifs perso
 383       hits[10] = mom[3]; // hit momentum P
 384       hits[11] = mom[0]; // Px/P
 385       hits[12] = mom[1]; // Py/P
 386       hits[13] = mom[2]; // Pz/P
 387       // fin modifs perso
 388
 389       // phi angle of incidence
 390       tlength = 0;
 391       eloss   = 0;
 392       eloss2  = 0;
 393       xhit    = pos[0];
 394       yhit    = pos[1];
 395       // Only if not trigger chamber
 396       if(idvol<10) {
 397           //
 398           //  Initialize hit position (cursor) in the segmentation model
 399           ((AliMUONChamber*) (*fChambers)[idvol])
 400               ->SigGenInit(pos[0], pos[1], pos[2]);
 401       } else {
 402           //geant3->Gpcxyz();
 403           //printf("In the Trigger Chamber #%d\n",idvol-9);
 404       }
 405   }
 406   eloss2+=destep;
 407
 408   //
 409   // Calculate the charge induced on a pad (disintegration) in case
 410   //
 411   // Mip left chamber ...
 412   if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
 413       gMC->SetMaxStep(big);
 414       eloss   += destep;
 415       tlength += step;
 416
 417       // Only if not trigger chamber
 418       if(idvol<10) {
 419           if (eloss > 0)
 420               MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),eloss,0.0,idvol);
 421       }
 422
 423       hits[6]=tlength;
 424       hits[7]=eloss2;
 425       if (fNPadHits > (Int_t)hits[8]) {
 426           hits[8]= hits[8]+1;
 427           hits[9]= (Float_t) fNPadHits;
 428       }
 429
 430       new(lhits[fNhits++])
 431           AliMUONHit(fIshunt,gAlice->GetCurrentTrackNumber(),vol,hits);
 432       eloss = 0;
 433       //
 434       // Check additional signal generation conditions
 435       // defined by the segmentation
 436       // model (boundary crossing conditions)
 437   } else if
 438       (((AliMUONChamber*) (*fChambers)[idvol])
 439        ->SigGenCond(pos[0], pos[1], pos[2]))
 440   {
 441       ((AliMUONChamber*) (*fChambers)[idvol])
 442           ->SigGenInit(pos[0], pos[1], pos[2]);
 443 //      printf("\n-> MakePadHits, reason special %d",ipart);
 444       if (eloss > 0)
 445           MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),eloss,0.0,idvol);
 446       xhit     = pos[0];
 447       yhit     = pos[1];
 448       eloss    = destep;
 449       tlength += step ;
 450       //
 451       // nothing special  happened, add up energy loss
 452   } else {
 453       eloss   += destep;
 454       tlength += step ;
 455   }
 456 }
 457
 458 //___________________________________________
 459 void AliMUONproto::BuildGeometry()
 460 {
 461     TNode *Node;
 462     TNode *Top;
 463
 464     const int kColorMUON = kBlue;
 465     //
 466     Top=gAlice->GetGeometry()->GetNode("alice");
 467     //
 468     //
 469     //
 470     float dx, dy, dz, zpos;
 471
 472     const Float_t cz[1]={975.};
 473
 474     zpos=cz[0];
 475
 476     dx=11*.975/2;
 477     dy=18*.55/2;
 478     dz=0.325;
 479
 480     new TBRIK("C_MUON101","Mother volume for Proto.","void",dx*2,dy*2,dz*2);
 481     TBRIK* PROTO = new TBRIK("PROT", "Proto. Sens. Region","void",dx,dy,dz);
 482     Top->cd();
 483     Node = new TNode("MUON101","ChamberNode","C_MUON101",0,0,zpos,"");
 484     Node->SetLineColor(kColorMUON);
 485     Node->SetVisibility(0);
 486     fNodes->Add(Node);
 487     Node->cd();
 488     Node = new TNode("MUON201", "Proto. Sens. Region Node", PROTO, -dx, dy, dz);
 489     Node->SetLineColor(kColorMUON);
 490
 491 }
 492
 493 //_____________________________________________________________________________
 494 void AliMUONproto::FindClusters(Int_t nev, Int_t last_entry)
 495 {
 496
 497 //
 498 // Loop on chambers and on cathode planes
 499 //
 500   for (Int_t icat = 0; icat < 2; icat++) {
 501             gAlice->ResetDigits();
 502             gAlice->TreeD()->GetEvent(last_entry+icat);
 503
 504       for (Int_t ich = 0; ich < kNTrackingCh; ich++) {
 505           AliMUONChamber* iChamber=(AliMUONChamber*) (*fChambers)[ich];
 506           TClonesArray *MUONdigits  = this->DigitsAddress(ich);
 507           if (MUONdigits == 0) continue;
 508           //
 509           // Get ready the current chamber stuff
 510           //
 511
 512           AliMUONResponse* response = iChamber->ResponseModel();
 513           AliMUONSegmentation* seg = iChamber->SegmentationModel(icat+1);
 514           AliMUONClusterFinder* rec = iChamber->ReconstructionModel();
 515
 516           if (seg) {
 517               rec->SetSegmentation(seg);
 518               rec->SetResponse(response);
 519               rec->SetDigits(MUONdigits);
 520               rec->SetChamber(ich);
 521               rec->FindRawClusters();
 522           }
 523
 524           TClonesArray *fRch;
 525           fRch=RawClustAddress(ich);
 526           fRch->Sort();
 527           // it seems to work
 528
 529       } // for ich
 530       // fill the tree
 531       gAlice->TreeR()->Fill();
 532
 533       ResetRawClusters();
 534   } // for icat
 535
 536   char hname[30];
 537   sprintf(hname,"TreeR%d",nev);
 538   gAlice->TreeR()->Write(hname);
 539   gAlice->TreeR()->Reset();
 540 }
 541
 542 #ifdef WE_FORGRET_THIS_SHIT
 543 //_____________________________________________________________________________
 544 void AliMUONproto::SetThreshold()
 545 {
 546
 547     ifstream inputFile("/home/alice/guernane/aliroot/pro/MUON/crped190.dat", ios::in);
 548
 549     if (inputFile.fail()) {
 550         AliError("Error opening file");
 551         exit(2);
 552     }
 553
 554     char buff[32];
 555     Int_t Serial;
 556     Int_t Ntrigger;
 557     Int_t Nchannel;
 558     Int_t i1;
 559     Int_t i2;
 560
 561     inputFile >> buff;
 562
 563     inputFile >> Serial;
 564     inputFile >> Ntrigger;
 565     inputFile >> Nchannel;
 566     inputFile >> i1;
 567     inputFile >> i2;
 568
 569     Float_t ped0[Nchannel];
 570     Float_t sig0[Nchannel];
 571     Float_t ped1[Nchannel];
 572     Float_t sig1[Nchannel];
 573     Int_t ichannel;
 574
 575     for (Int_t i = 0; i < Nchannel-1; i++) {
 576     ped0[i] = 0;
 577     sig0[i] = 0;
 578     ped1[i] = 0;
 579     sig1[i] = 0;
 580     }
 581
 582     for (Int_t i = 0; i < Nchannel-1; i++) {
 583         inputFile >> ichannel;
 584         inputFile >> ped0[i];
 585         inputFile >> sig0[i];
 586         inputFile >> ped1[i];
 587         inputFile >> sig1[i];
 588         fThreshold[i] = fNsigma*sig1[i];
 589     }
 590
 591     inputFile.close();
 592 }
 593 #endif